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ABSTRACT 

Nonlinear evolution causes the galaxy power spectrum to become broadly corre- 
lated over different wavenumbers. ft is shown that prewhitening the power spectrum 
- transforming the power spectrum in such a way that the noise covariance becomes 
proportional to the unit matrix - greatly narrows the covariance of power. The eigen- 
functions of the covariance of the prewhitened nonlinear power spectrum provide a set 
of almost uncorrelated nonlinear modes somewhat analogous to the Fourier modes of 
the power spectrum itself in the linear, Gaussian regime. These almost uncorrelated 
modes make it possible to construct a near minimum variance estimator and Fisher 
matrix of the prewhitened nonlinear power spectrum analogous to the Feldman-Kaiser- 
Peacock estimator of the linear power spectrum. The paper concludes with summary 
recipes, in gourmet, fine, and fastfood versions, of how to measure the prewhitened 
nonlinear power spectrum from a galaxy survey in the FKP approximation. An Ap- 
pendix presents FFTLog, a code for taking the fast Fourier or Hankel transform of a 
periodic sequence of logarithmically spaced points, which proves useful in some of the 
manipulations. 
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1 INTRODUCTION 

Most of the information about cosmological parameters bot- 
tled inside current^ and coming galaxy surveys, notably the 
Two-Degree Field Survey (2dF) (Colless 1998; Folkes et al. 
1999) and the Sloan Digital Sky Survey (SDSS) (Gunn & 
Weinberg 1995; Margon 1998), lies in the nonlinear regime. 
Even in the linear regime, nonlinearities perturb. 

At large, linear scales, the power spectrum - the covari- 
ance of the density field, expressed in the Fourier representa- 
tion - is the preeminent measure of large scale structure. It 



* For a review of redshift surveys of galaxies see Strauss (1999) 
and references therein. Recent surveys include: Updated Zwicky 
Catalog (UZC) (Falco et al. 1999); IRAS Point Source Catalogue 
Redshift Survey (PSCz) (Sutherland et al. 1999); Redshift Survey 
of Zwicky Catalog Galaxies in a 2'' X 15 deg Region around 3C 273 
(Grogin, Geller & Huchra 1998); Durham/UKST (Ratcliffe et al. 
1998); Southern Sky Redshift Survey (SSRS) (da Costa et al. 
1998); ESQ Slice Project (ESP) (Vettolani et al. 1998); Muen- 
ster Redshift Project (MRSP) (Schuecker et al. 1998); CNOC2 
Field Galaxy Redshift Survey (Carlberg et al. 1998); Century 
Survey (Geller et al. 1997); Norris Survey of the Corona Borealis 
Supercluster (Small, Sargent & Hamilton 1997); Stromlo-APM 
(Loveday et al. 1996); Las Campanas Redshift Survey (LCRS) 
(Shectman et al. 1996); Hawaii Deep Fields (Cowie et al. 1996); 
Canada- France Redshift Survey (CFRS) (Lilly et al. 1995). 



is a generic, though by no means universal, prediction of in- 
flation (Turner 1997) that linear density fluctuations should 
be Gaussian. More generally, primordial fluctuations should 
be Gaussian whenever they result from superpositions of 
many independent processes, thanks to the central limit the- 
orem. Observations of large scale structure are consistent 
with linear density fluctuations being Gaussian (Bouchet et 
al. 1993; Juszkiewicz, Bouchet, & Colombi 1993; Gaztanaga 
1994; Gaztanaga & Frieman 1994; Nusser, Dekel & Yahil 
1995; Stirling & Peacock 1996; CoUey 1997; Chiu, Ostriker 
& Strauss 1998; Frieman & Gaztanaga 1999) although the 
evidence is not deflnitive (White 1999). If linear density 
fluctuations are Gaussian, then the 3-point and higher irre- 
ducible moments are zero, so that the covariance of the den- 
sity field contains complete information about the statistical 
properties of the field, hence all information about cosmo- 
logical parameters. Compared to other measures of covari- 
ance such as the correlation function, the power spectrum 
has the additional advantage that estimates of power at dif- 
ferent wavenumbers are uncorrelated, for Gaussian fluctua- 
tions. This asset of the power spectrum is intimately related 
to the assumption that the field is statistically translation 
invariant, and to the fact that Fourier modes are eigenfunc- 
tions of the translation operator. 

At smaller, nonlinear scales, the power spectrum loses 
some of its glow. Nonlinear evolution drives the density 
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field away from Gaussianity, coupling Fourier modes, feed- 
ing higher order moments, and causing power at different 
wavenumbers to become correlated. The broad extent of the 
correlation of the nonlinear power spectrum has been em- 
phasized by Meiksin & White (1999) and Scoccimarro, Zal- 
darriaga & Hui (1999), and is illustrated in Figure ^of the 
present paper. 

The purpose of the present paper is to show how to un- 
fold the nonlinear power spectrum into a set of nearly uncor- 
related modes, somewhat analogous to the Fourier modes of 
the power spectrum itself in the linear, Gaussian regime. The 
present paper is a natural successor to Hamilton (1997a,b, 
hereafter Papers 1 and 2), which showed how to derive the 
minimum variance estimator and Fisher matrix of the power 
spectrum of a galaxy survey in the Feldman, Kaiser & Pea- 
cock (1994, hereafter FKP) approximation, for Gaussian 
fluctuations. Section 5.2 of Paper 1 posed, but was unable 
to solve, the non-Gaussian problem solved in the present 
paper. A following paper (Hamilton & Tegmark 2000, here- 
after Paper 4), describes how to complete the processing of 
the power spectrum into fully decorrelated band-powers. 

It turns out that a key to solving the non-Gaussian 
problem is to 'prewhiten' the power spectrum - to trans- 
form the nonlinear power spectrum in such a way that the 
(2-point) shot-noise contribution to the covariance matrix is 
proportional to the unit matrix. The properties of the pre- 
whitened nonlinear power spectrum appear empirically to 
be sweeter than might reasonably have been expected. 

This paper is devoted entirely to the problem of nonlin- 
earity. It ignores the equally important problem of redshift 
distortions (Hamilton 1998), and the problematic question 
of light-to-mass bias (Coles 1993; Fry & Gaztanaga 1993; 
Mo, Jing & White 1997; Mann, Peacock & Heavens 1998; 
Tegmark & Peebles 1998; Moscardini et al. 1998; Scherrer 
& Weinberg 1998; Dekel & Lahav 1999; Colin et al. 1999; 
Gen & Ostriker 1999; Narayanan, BerUnd & Weinberg 1999; 
Blanton et al. 1999; Benson et al. 1999; Bernardeau & Scha- 
effer 1999; Coles, Melott & Munshi 1999). It further assumes 
that uncertainties arising either from the selection function 
(Binggeli, Sandage & Tammann 1988; Willmer 1997; Tresse 
1999) or from evolution in the cosmological volume element 
or the gala:xy population, are negligible. 

Several authors have recently published estimates of 
how well measurements of the power spectrum from fu- 
ture galaxy surveys will constrain cosmological parameters 
(Tegmark 1997b; Goldberg & Strauss 1998; Hu, Eisenstein 
& Tegmark 1998; Eisenstein, Hu & Tegmark 1998, 1999). 
The procedures described in the present paper should assist 
this enterprise. 

The aims of the present paper are complementary to 
those of Bond, Jaff'e & Knox (1998b). The question Bond 
et al. considered was: If the power spectrum (of the Cosmic 
Microwave Background, specifically) is quadratically com- 
pressed (Tegmark 1997a; Tegmark et al. 1997, 1998) into a 
set of band-powers, then what is the best way to use those 
band-powers in Maximum Likelihood estimation of param- 
eters? For example, one general procedure is to use not the 
band-powers themselves, but rather functions of the band- 
powers arranged such that their variances remain constant 
as the prior power is varied. Bond et al. argued that the 
likelihood function is then more nearly Gaussian. The pur- 
pose of this paper and Paper 4 is rather to arrive at the 



point where one has decorrelated band-powers to work with 
in the first place. 

The plan of this paper is as follows. Section ^ sets up the 
notation and defines reference material needed in subsequent 
sections. Section ^ goes through the difficulties one meets 
in attempting to measure the nonlinear power spectrum in 
minimum variance fashion, and describes how to overcome 
them. Section ^ reveals the unexpectedly nice properties of 
the prewhitened covariance of the power spectrum, key to 
the whole enterprise of this paper. Section |5 defines the 
prewhitened power spectrum. Sections ^ and ^ show how 
the approximations motivated in previous sections lead to 
a practical way to evaluate the Fisher matrix of the pre- 
whitened nonlinear power, and to measure the prewhitened 
nonlinear power spectrum from a galaxy survey. Section ^ 
discusses how to evaluate the Fisher matrix and nonlinear 
power spectrum using the FKP approximation alone, with- 
out any additional approximation. Section |^ summarizes the 
results of previous sections into recipes, in gourmet, fine, and 
fastfood versions, for measuring nonlinear power, the end 
product being a set of uncorrelated prewhitened nonlinear 
band-powers, with error bars, over some prescribed grid of 
wavenumbers. Section |l^ summarizes the conclusions. Ap- 
pendix B gives details of FFTLog, a code for taking the 
fast Fourier or Hankel transform of a periodic sequence of 
logarithmically spaced points. 



2 PRELIMINARIES 

This section contains reference material needed in subse- 
quent sections. The reader interested in new results may 
like to skip to the next section, f|^, referring back to the 
present section as needed. 

2.1 Data, parameters 

'He will, of course, use maximum likelihood because his text- 
books have told him that' - E. T. Jaynes (1996, p. 624). 

According to Bayes' theorem, the probability distribu- 
tion of parameters 6c, given data Si is, up to a normalization 
factor, the product of the prior probability with the likeli- 
hood function £.{Si\9a). The data Si in a galaxy survey can 
be taken to be overdensities S{r) at positions r in the survey 

n{r) — n{r) 



n(r) 



(1) 



where n(r) is the observed number density of galaxies, and 
n(r) is the selection function. The parameters Oa are, for the 
present purpose, some parametrization of the galaxy power 
spectrum; the focus of this paper is on the case where the 
parameters are the power spectrum ^c, itself. 

This paper conforms to the common convention used 
by cosmologists to relate the power spectrum ^(fc) in Fourier 
space to the correlation function ^(r) in real space, notwith- 
standing the extraneous factors of 27r that result: 



''^{r)dr— / jo{kr)£,{r) Atvt dr 
Jo 



io(fcr)C(fc) 



47rfc dfc 



(2^)3 ' (2^)3 

where jo{x) — sinx/x is a spherical Bessel function. 



(2) 
(3) 
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2.2 Hilbert space 

As in Paper 1, it is convenient to adopt a notation in which 
Latin indices i, j, refer to 3-dimensional positions, while 
Greek indices a, (3, run over the space of parameters, and 
more specifically over the 1-dimensional space of wavenum- 
bers or pair separations. 

For generality, brevity, and ease of manipulation, it is 
convenient to treat quantities such as the data vector Si , or 
the power spectrum ^a, as vectors in a Hilbert space (for 
a didactic exposition, see Hamilton 1998 §3.3). Such vec- 
tors have a meaning independent of the particular basis, i.e. 
complete set of linearly independent functions, with respect 
to which they might be expressed. For example, the data 
vector has components 5r [— S{r)] when expressed in real 
space, or components Sk [= S{k) = J e^'°'^S{r) dV] when ex- 
pressed in Fourier space, but from a Hilbert space point of 
view these are the same vector, and in this paper they are 
both denoted by the same symbol Si. 

Similarly the power spectrum has components 
[— 5(fc)] when expressed in Fourier space, or [~ £,{r)] when 
expressed in real space, but again from a Hilbert space point 
of view these are the same vector, and in this paper they are 
both denoted by the same symbol ^a- 

Latin indices i, j, on vectors and matrices run over 
the 3-dimensional space of positions r, or more generally 
over any 3-dimensional basis of the Hilbert space. Unless 
stated otherwise, repeated pairs of indices signify the inner 
product in Hilbert space, as in 



abi 



a*{r)b{r) d r - 



a*{k)b{k) d\/(2nf 



(4) 



By definition, the inner product is a scalar, the same quan- 
tity independent of the choice of basis. The raised index a* 
denotes the Hermitian conjugate (if the basis is orthonor- 
mal) of the vector Ui. One of the indices in an inner product 
is always raised, the other lowered. In this paper, all vectors 
in the Hilbert space are real-valued when expressed in real 
space, so that a*(r) = a{r) and a*(fe) — a{—k). 

Adhering to the raised/lowered index convention serves 
as a useful reminder that one of the pair of vectors in an 
inner product is a Hermitian conjugate (if the basis is or- 
thonormal). In Fourier space, for example, this means using 
—k for one index (raised) and +k for the other index (low- 
ered) of an inner product. 

Greek indices a, (5, run over the space of 1- 
dimensional pair separations r, or wavenumbers k, or more 
generally over any 1-dimensional basis in the associated 
Hilbert space. Again, unless stated otherwise, repeated in- 
dices signify the inner product 



a*{r)b(r)4nr^dr ^ / a*{k)b{k) 4nk^dk/{2Tvf (5) 



a ba 



which is again a scalar, the same quantity independent of 
the choice of basis. Again, in this paper all vectors in the 
Hilbert space are real-valued in real space, so a*{r) — a{r) 
and a*{k) = a{k). Although there is no distinction in this 
case between vectors with raised and lowered indices in ei- 
ther real or Fourier space, adhering to the raised/lowered 
index convention again serves as a useful reminder. 

The unit matrix if in any representation is defined such 



that its inner product with any vector ap leaves the vector 
unchanged, 

la CL0 = ap I'l = Ua ■ (6) 

In the continuous real representation, the unit matrix is 

Irf = S'ioira - f/s) (7) 

where S-inira — rp) denotes the 3-dimensional Dirac delta- 
function, defined such that 



S'ioira - rp) 47rr^dra = 1 



(8) 



In the continuous Fourier representation, the unit matrix is 
lll^i27vfS3D{ka^kf>) (9) 
again a 3-dimensional Dirac delta-function. 



2.3 Discretization of matrices 

Many of the operations in this paper involve manipulations 
of matrices in the 1-dimensional space of separations. Con- 
tinuous matrices must be discretized to manipulate them 
numerically. Discretization should be done in such a way as 
to preserve the inner product (^), so that integration over the 
volume element, Anr^dr in real space, or 47vk'^dk/{2TTf m 
Fourier space, translates into summation in the correspond- 
ing discrete space. This ensures that matrix operations such 
multiplication, diagonalization, and inversion can be done 
in the usual fashion. 

Most of the manipulations in this paper are done in 
Fourier space on a logarithmically spaced grid of wavenum- 
bers ka. In this case, a continuous vector a{kc,) is discretized 
by multiplying it by [47rfe^ A In fc/(27r)^]^''^ 

a(fc^) afe„ = a(fc„) [47rfci A In k/{2T^f] '''^ (10) 

and a continuous matrix A{ka,kfj) is discretized by multi- 
plying it by 47r(fc„fc/3)^''^Alnfc/(27r)^ 

A{kc,,k0) ^ kk^k^ = A{kc,,ki3)4n{kak0f/^Alnk/{27Tf .(11) 

The unit matrix (27r)'^53_D(fca — fc^j) in the continuous Fourier 
representation translates to the unit matrix 1^/3 in the dis- 
crete case 



i27vYS3D{ka-kp) 



^13 



(12) 



Similarly, a continuous vector aire) in real space is dis- 
cretized on to a logarithmically spaced grid of separations 
Ta by multiplying the vector by [47rr|^Alnr] ^^'^ 

a{ra) —* a,.„ = a{ra) [47rr-^ A In r] "'"''^ (13) 

and a continuous matrix A{ra,r0) is discretized by multi- 
plying it by 4Tv{rarp)^^^ Alnr 

A(ra,r^) Ar^r^ = A{ra,r0)4Tv{rc,rf3f/'^Alnr . (14) 

The unit matrix Saoirci — rp) in the continuous real repre- 
sentation translates to the unit matrix la/3 in the discrete 
case 



S'ioira - r(3) 



(15) 



The transformation between Fourier and real space for 
logarithmically spaced wavenumbers ka and separations ra 
may be accomplished with FFTLog (Appendix B). 
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2.4 Gaussian density field 



(26) 



If the density distribution S{r) were Gaussian - wliicli is not 
true in tlie present case - tiien one would liave tiie luxury 
of being able to write down an explicit Gaussian likelihood 
function 

1 .expf-i5.C-"^5,) (16) 



C oc 



|C|l/2 , _ 

where \C\ and are the determinant and inverse of the 
covariance matrix C of overdensities 



(17) 



Angle-brackets here and throughout this paper signify av- 
erages over possible data sets Si predicted by the likelihood 
function 



(18) 



Maximum Likelihood (ML) estimates 9a of the parameters 
9a (the hat distinguishing the estimate 9a from the true 
value 9a) are given by the vanishing of the vector of partial 
derivatives of the log-likelihood function 

dlnC Ida- 



d9a 
dlnC 



2 961, 



'-C- 



d9a 



(19) 
(20) 



The covariance {A9aA9i3) of the estimated parameters is 
given approximately by the inverse of the Fisher information 
matrix defined to be minus the expectation value of 

the matrix of second partial derivatives of the log-likelihood 
function 



F 



al3 



d9ad9f, 



{A9aA9p) ^ Fj . 



1 dCij 



d9g 



(21) 



(22) 

The approximation (^^ is exact if the estimated parameters 
9a are Gaussianly distributed about their expectation val- 
ues. The central limit theorem asserts that the parameters 
become Gaussianly distributed in the asymptotic limit of a 
large amount of data. 

It is commonly assumed, and the same assumption is 
adopted here, that the dominant source of variance in a 
galaxy survey is a combination of cosmic (sample) variance 
and shot-noise arising from the discrete sampling of galax- 
ies. If the sampling of galaxies is random - a Poisson process 
- then the covariance dj is a sum of the cosmic covariance 



^ij with Poisson sampling noise Ni- 
Cij = -|- Nij . 



(23) 



In the real representation, the cosmic covariance is the 
correlation function 



and the noise matrix Nij is the diagonal matrix 



N.J = 



5'iD{ri - Tj) 
n{ri) 



(24) 



(25) 



with 5iD{ri — Vj) a 3-dimensional Dirac delta- function. In 
the Fourier representation the cosmic covariance is the 
diagonal matrix 



whose eigenvalues €,{ki) constitute the power spectrum. 

The focus of this paper is on the case where the param- 
eters 9a are the power spectrum itself (in this paper the 
cosmic covariance function , expressed in an arbitrary rep- 
resentation, will often be referred to as the 'power spectrum', 
even though this name is commonly reserved for the covari- 
ance ^(fc) expressed in Fourier space; no confusion should 
result). In this ca,s6 the covariance Cij is a linear function of 
the parameters 



Cij — Dij^a ~\~ ^ij 



(27) 



where in real space £,a = 5('"a) is the correlation function, 
and 

DT,^S3D{\r,-rj\-ra) (28) 

is a 3-dimensional Dirac delta- function, equation (|l^), while 
in Fourier space ^a = C(^q) is the thing commonly called 
the power spectrum, and 



(29) 



It follows from equations ( p^ and (^o|) that the ML es- 
timator ^Q, of the power spectrum, for Gaussian fluctuations, 
is that solution of 



ic 



for which the estimate is equal to the prior, 
variance of the ML estimator is 

{AiaAip) ^ F-^ 

and the Fisher matrix is 



(30) 
^a- The 



(31) 



(32) 



If the prior power is regarded as fixed, then equa- 
tion (|^) yields an estimated power that is quadratic in 
overdensities Si. If this estimated power is folded back into 
the prior, then equation (^) with the revised prior yields 
another estimate of power. Iterated to convergence, the re- 
sult is the ML estimator of the power. It is to be noted that 
even without iteration, equation (^^ yields a measurement 
of power that (as long as the prior is at least roughly correct) 
should already be a good approximation, since 'if the prior 
matters, then you are not learning much from the data', to 
quote one of the refrains from the 1997 Aspen workshop on 
Precision Measurement of Large Scale Structure. 

The question of how to apply quadratic estimators (such 
as given by equation [^) to measure the power spectrum 
is addressed by Tegmark et al. (1998) for galaxies, and by 
Tegmark (1997a), Tegmark et al. (1997), and Bond, Jaffe & 
Knox (1998a,b) for the CMB. 



2.5 Non-Gaussian density field 

Ultimately, one might look forward to a wondrous A'^-body 
machine able to compute the probability distribution of lin- 
ear initial conditions given noisy and incomplete data from a 
survey (Narayanan & Weinberg 1998; Monaco & Efstathiou 
1999; and references therein). 

In the meantime it is far from clear what to write down 
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as a likelihood function for the nonlinear density field (Do- 
delson, Hui & Jaffe 1999). Certainly it would be a bad idea 
to use a Gaussian likelihood function for a non-Gaussian 
density field, since that would lead to a serious underes- 
timate of the true uncertainty in the measured nonlinear 
power spectrum. 

An alternative procedure is to seek a minimum variance 
unbiased estimator of power. Now the power spectrum is 
by definition a covariance of overdensities, and by the pre- 
sumption of Poisson sampling, any a priori weighted sum 
of quantities quadratic in observed overdensities (with self- 
terms excluded, to eliminate shot-noise) provides an unbi- 
ased estimate of the power spectrum linearly windowed in 
some fashion. It was shown in §2.3 of Paper 1 that, amongst 
estimators quadratic in observed overdensities Si, the unbi- 
ased estimator ^„ of the power spectrum having minimum 
variance is 



- ap ij 

with variance 

where F"^ is the Fisher matrix 

_ n« — 'i-ijkl p.f3 



(33) 



(34) 



(35) 



^ijki is the covariance of shot-noise-subtracted products of 
overdensities 



{{SiSj - Nij - Cij) (SkSi ~ Nki - CfcO) 



(36) 



and is its inverse, meaning CijmnC"^™"*' = 

Sym^^jjjl^lj. The symbol Sym^j^j signifies symmetrization 
over its underscripts, as in 



Sym Aij = ^{A^J+Aji) 

(ij) ^ 



(37) 



The quantity Nki in equations (^3|) and ( ^ ) is the 'ac- 
tual' shot-noise, the contribution to S^Si from self-pairs of 
galaxies, pairs consisting of a galaxy and itself. The actual 
shot-noise Nm in a survey is to be distinguished from its 
expectation value A^^; = (Nki)- If the expected shot-noise 
Nki is used in equation (KSl) in place of the actual shot- 
noise, then additional terms (given in eq. [8] of Paper 1) 
appear in the covariance matrix <Zijki, increasing the vari- 
ance of the estimator. Why does the ML estimator in 
the Gaussian case, equation (^) , involve the expected shot- 
noise Nki rather than the actual shot-noise Afei? Because 
a discretely sampled Gaussian field is not really Gaussian, 
except in the limit where a cubic wavelength contains many 
galaxies, so the assumption of a Gaussian likelihood function 
is not strictly correct. In fact it is plain that the Gaussian 
ML estimator would also be improved if the actual shot- 
noise Nki were used in place of the expected shot-noise Nki 
in equation (|3C|), since using the actual shot-noise exploits 
additional information about the character of the Poisson 
sampling that is discarded by the Gaussian likelihood. How- 
ever, as discussed by Tegmark et al. (1998 Appendix A), 
the gain from subtracting the actual versus the expected 
shot-noise is in practice small at linear scales, where a cubic 
wavelength is likely to contain many galaxies. 

In the same Poisson sampling approximation as equa- 
tion (Ei), the covariance ^ijki of sliot-noise-subtracted prod- 





Figure 1. Schematic illustration of the 4-point, 3-point, and 2- 
point contributions to the covariance €.ijkl of pairs ij with other 
pairs kl. The 3-point and 2-point contributions arc shot-noise 
contributions in which one or both galaxies of the pair ij are the 
same as one or both of the pair kl. 



nets of overdensities, equation (|3q), is, in the real represen- 
tation with no implicit summation. 



+ [Niki^ji + Oji) + ii^ j,k 
+ {N.kNji + NuNjk){l + C»j) 



;)] (4 terms) 



(38) 



in which the top line is the 4-point, the middle the 3-point, 
and the bottom line the 2-point contribution to the covari- 
ance, as illustrated in Figure ^. For Gaussian density fiuc- 
tuations equation (^) reduces to 

(i,,ki = 2 Sym CkCji (39) 
(kl) 

with inverse 

^f-Ujkl ^ 1 ^- 

2 (kl) 

It follows from equation ( ]4o| ) that for Gaussian fluctua- 
tions the minimum variance estimator of the power spec- 
trum, equation (^), is the same as the ML estimator, equa- 
tion (Isd), if the estimate is folded back into the prior and 
iterated to convergence (modulo the comments about shot- 
noise in the previous paragraph). 



(40) 



3 PROBLEMS 

3.1 FKP approximation 

Calculating the minimum variance estimate of the power 
spectrum, equation (^), involves the formidable problem of 
inverting the pair covariance ^ijki, a rank 4 matrix of 3- 
dimensional quantities. Whereas for Gaussian fluctuations 
the rank 4 matrix <Liiki factorizes into a product of rank 
2 matrices, equation (|39[), for non- Gaussian fluctuations it 
does not factorize. Again, whereas for Gaussian fluctuations 
it may be possible, at least at the largest scales, to pixelize 
a survey into large enough pixels that brute force numeri- 
cal inversion is feasible, for non-Gaussian fluctuations brute 
force inversion is quite impossible. 

A natural way to simplify the problem is to adopt the 
Feldman, Kaiser & Peacock (1994, FKP) approximation, 
where the selection function n(r) is taken to be locally con- 
stant. The FKP approximation is expected to be valid at 
wavelengths much smaller than the characteristic size of the 
survey. Section 5 of Paper 1 terms this the 'classical' ap- 
proximation, since it is valid to the extent that the position 
and wavelength of a density mode can be measured simul- 
taneously. While the FKP approximation is liable to break 
down at larger scales, particularly for pencil beam or slice 
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surveys, it should be a good approximation at smaller, non- 
linear scales, especially in surveys with broad contiguous sky 
coverage. 

Even if the selection function n is taken to be constant, 
the general problem of inverting the rank 4 matrix (Lijki 
remains intractible. Notice however that g^"!*^*' appears 
multiplied in both equations (^) and ( p5[ ) by the matrix 
Dfj. Now Dfj has translation and rotation symmetry, and 
in the FKP approximation the matrix Cijti also has transla- 
tion and rotation symmetry, the selection function n being 
constant. Indeed, inspection of equation reveals that 
the matrix ^ijki remains translation and rotation invariant 
even if the selection functions and fij at positions i and 
j are two different constants. It follows that the combina- 
tion >ZijkiD^ is likewise translation and rotation symmetric, 
which implies that it can be expressed in the form 



€ai3{ni,nj) D. 



(41) 



for some matrix £a/3, which can be termed the 'reduced' 
covariance matrix. Equation (^]) is the FKP approxima- 
tion, expressed in concise mathematical form; additional de- 
tails of the justification of this equation are provided in Ap- 
pendix A. The reduced matrix is written in equation (^l|) as 
<tap{ni,nj) to emphasize the fact that it is a function of the 
selection functions fii and fij at positions i and j; note that 
no implicit summation over i or j is intended on the right 
hand side of equation (^). Inspection of equation ( psj ) for 
Cijfci shows that the reduced covariance <Lap{ni,nj) takes 
the form 

<Zapini,nj) = 2 [KcP + {n~^+nJ^)Jap + n~^n~^ Hap] (42) 

a linear combination of 4-point, 3-point, and 2-point contri- 
butions Kap, Jap, and Hap. Multiplying equation ( ^l| ) by 
ff'-i-ya ^f-imnij gjjg^g ^j^g^^ ^^le iuverse of €ijki is similarly 
related to the inverse of the reduced matrix <tap 



-lijkl j-^a 



-lap,- - 



{rii, 71^)0'^ 



(43) 



Physically, to the extent that the selection functions fii 
and fij at positions i and j are constants, the minimum vari- 
ance pair-weighting attached to a pair ij should be a func- 
tion only of the separation a of the pair, not of their position 
or orientiation. Just as Cijfc; is the covariance between a pair 
ij and another pair kl, so the reduced covariance matrix C^/s 
is the covariance between a pair separated by a and another 
pair separated by /3. 

In the FKP approximation given by equation (^), the 
minimum- variance estimate ( ^ ) of the power spectrum is 

ia = F-;€-^^-'Dl^\5,5, - N,,) (44) 



and the associated Fisher matrix (1351) is 



F 



ap _ — la7 rfi n/3 



(45) 



Notice that the approximate Fisher matrix given by this 
equation (^) is not symmetric, whereas the original Fisher 
matrix, equation (psj), was symmetric. The asymmetry re- 
sults from the asymmetry of the FKP approximation, equa- 
tion (^). The approximate expression ( |45[ ) would be sym- 
metric if the FKP approximation were exact, and in practice 
it should be nearly symmetric; if not, it is a signal that the 
FKP approximation is breaking down. 

To ensure symmetry of the Fisher matrix, one might 
be inclined at this point to symmetrize equation (^), since 



after all an equally good approximation to the Fisher matrix 
would be the same expression ( p5| ) with the indices swapped 
on the right hand side, a « /3. However, it is desirable that 
the FKP estimator equation (^), should be unbiased, 
meaning that 



Averaging equation 
according to equation 



gives, since 



(Co) 



^ap^ -^7 



(46) 
(47) 



which shows that the FKP estimator is unbiased only if 



is interpreted as satisfy- 
A detailed discussion of 



the Fisher matrix in equation i 
ing the asymmetric expression 
this issue is deferred to Here it suffices to remark that, 
to the extent that the FKP approximation is valid, the vari- 
ance of the FKP estimator is equal to the inverse of the 
symmetrized Fisher matrix given by equation (143) 



where F'^"' 



(aP) 



(48) 



Sym,^0^F°''^ denotes the symmetrized Fisher 



matrix, and F, \, its inverse. 

' (a/3) 



3.2 Hierarchical model 



hence also 
4l|), involves 



The pair covariance matrix ^ijki, equation (|38D 
the reduced covariance matrix <Lap, equation 
the 3-point and 4-point correlation functions Qjk and rjijki- 
The problem here is that these correlation functions are not 
known precisely. 

Available observational and A''-body evidence (see for 
example the summaries by Scoccimarro & Frieman 1999 
and Hui & Gaztanaga 1999) is consistent with a hierarchical 
model in which the 3-point and 4-point functions are, in the 
real representation with no implicit summation. 



(49) 



Vijki = i?a [CijCjfcCfci + cyclic (12 snake terms)] 

+ Rb [^ij^ik^ii + cyclic (4 star terms)] (50) 

with approximately constant hierarchical amplitudes Q, Ra, 
and Rb- On the other hand it is clear that the hierarchical 
amplitudes do vary at some level, both as a function of scale 
and configuration shape. 

In the translinear regime, perturbation theory predicts 
that the hierarchical amplitudes should vary (somewhat) 
with both scale and configuration, for density fluctuations 
growing by gravity from Gaussian initial conditions (Fry 
1984; Scoccimarro et al. 1998). 

In the deeply nonlinear regime, predictions for the be- 
haviour of the hierarchical amplitudes are more empirical. 
Scoccimarro & Frieman (1999) have recently suggested an 
ansatz, which they dub hyperextended perturbation theory 
(HEPT), that the hierarchical amplitudes in the highly non- 
linear regime go over to the values predicted by perturba- 
tion theory for configurations coUinear in Fourier space. For 
power law power spectra ^(fc) oc k", HEPT predicts a 3- 
point amplitude 



= 



- 2 2" 



(51) 
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and 4-point amplitudes Ra — Rb = Qi with 
54 - 27 2" + 2 3" + 6" 



(-Ra + Rb) 



Qa = 



2 (1 + 6 2" +3 3" + I 



(52) 



For simplicity, the present paper adopts the hierarchi- 
cal model, with constant hierarchical amplitudes set equal to 
the HEPT values (jHl]) and (|5^). For reasons to be discussed 
shortly (namely that the Schwarz inequality is violated oth- 
erwise), most of the calculations shown take 

Ra = --Rb = Qi (53) 

although where possible results are also shown for 



Ra — Rb — ' 



(54) 



In addition to power law power spectra, the present paper 
shows results for the power spectrum derived from observa- 
tions by Peacock (1997), and for an observationally concor- 
dant ACDM model from the fitting formulae of Eisenstein & 
Hu (1998), nonlinearly evolved according to the procedure of 
Peacock & Dodds (1996). In these cases the adopted ampli- 
tudes are those corresponding to n = —1.2, i.e. a correlation 
function with slope 7 — n + 3 — 1.8, for which Q = 1.906 
and Q4 = 4.195. 

In the hierarchical model with constant hierarchical am- 
plitudes, the 4-point, 3-point, and 2-point contributions to 
the reduced covariance matrix £a/3, equation (^), are, in 
the Fourier representation with no implicit summation, 

+ Ra[€.{k^)+£.{kfi)fA{k^,kfi) 
+ Rbakc)ak0mkc)+^{kp)] 



(55) 



J{ka,kf}) = (27r) S-ioika-kfj) £,{ka) 

+ QK(fcc.) + e(fc/3)]^(fcc,fc^) 

+ QC(fcc)C(fc/5) 

I{ka, kis) = {2l\)^53D{ka~kp) + A{ka,kf3) 



(56) 
(57) 



where in the real space representation the matrix Aa/s is the 
diagonal matrix 



A{ra,ri3) = Ssoira-r/s) ^(ra) 

while in the Fourier representation Aa/s is 



A{ka,ki3) 



2 k ^ 



£,(k) kdk 



(58) 



(59) 



Convergence of A{ka,kfs) at ka = kp requires that 
5(fc) ~ fc" with n > —2 at small wavenumber k. Conver- 
gence of J ^{r) Airr^dr at small r requires that ^(r) ~ r'"' 
with 7 < 3 at small separation r. Thus for power law power 
spectra ^{k) oc fc" (this is the evolved, nonlinear power spec- 
trum, not the original, linear power spectrum), equivalent to 
power law correlation functions ^(r) oc r ' with 7 = n -I- 3, 
the hierarchical model is consistent only for 



2 < n < or equivalently 1 < 7 < 3 



(60) 



It is straightforward to determine that, for power law 
power spectra £,{k) oc k" in the hierarchical limit (where the 
Gaussian contribution becomes negligible), the correlation 
coefficient of the 4-point contribution to the reduced 
covariance Cq/j is, for ka ^ kp, 



{KaaKfipy/^ 



2'i+2 
2 I ~nRa + Rb 



(fee) 



n/2 



(61) 



which diverges as ka/kp —> 00 (for —2 < n < 0) unless 
Rb — —Ra- Thus the Schwarz inequality, which requires that 
the absolute value of the correlation coefficient be less than 
or equal to unity, is violated unless Rb = — -Ra. This problem 
has been remarked and discussed by Scoccimarro et al. (1999 
§3.3). Scoccimarro et al. show from -/V-body simulations that 
the traditional relation -Ra « Rb holds approximately for 
kc, ~ kp, but that indeed Ra+Rb decreases systematically as 
ka and fc^ become more and more separated. Scoccimarro et 
al. conclude that the simple hierarchical model with constant 
amplitudes is not a good description of the 4-point function 
in the highly nonlinear regime. 

For simplicity, the present paper adopts the hierarchical 
model with constant amplitudes, and either Rb — —Ra or 
Rb = Ra- Ultimately, the latter choice leads to unphysically 
huge variances, plainly a consequence of the violation of the 
Schwarz inequality. Thus the canonical models in this pa- 
per have Rb = —Ra- However, where possible, intermediate 
results are also shown for -R;, — Ra- 



3.3 Prewhitening 

The minimum variance estimator f a and associated Fisher 
matrix F°'^ , equations (^) and (^5|), involve 6-dimensional 
integrals of ^-^"'^(n, , Uj ) over all pairs ij of volume ele- 
ments in a survey. This is actually quite a feasible numeri- 
cal problem. The reduced covariance matrix (tai3{ni,nj) is a 
rank 2 matrix of 1-dimensional quantities, so is straightfor- 
ward to invert numerically for any particular values of the 
selection functions fii and fij. If, as is typical, the selection 
function separates into the product of an angular mask and 
a radial selection function, then the angular integrals can be 
done analytically (Hamilton 1993), leaving a double integral 
of £-i"'3(n„n,) over the radial directions, which is doable. 
This direct procedure is discussed further in jjgl and forms 
the basis of the gourmet recipe summarized in j: 9jJ. Still, the 
integration is burdensome, and it is enlightening to explore 
whether further simplification is possible. 

Ideally what one would like is that there would exist a 
representation in which ltap{ni,nj) were simultaneously di- 
agonal for arbitrary values of the selection function n. Pre- 
cisely this situation obtains in the case of Gaussian fluc- 
tuations, for which the reduced covariance matrix (Lap is 
diagonal in Fourier space 

<^ap{n^,nj) =2{27vfS3D{ka-ki3)[^{ka)+fi-^] (fcc< ) (62) 

regardless of the values fii and fij of the selection function. 

For non-Gaussian fluctuations, the reduced covariance 
'tapifiijfij) is a linear combination of 4-point, 3-point, and 
2-point matrices Ka/s, Jap, and Ha/s, according to equa- 
tion (ji^). Finding a representation in which £^^(71,;,%) is 
diagonal for any fii and fij, thus means diagonalizing the 
three matrices K, J, and H simultaneously. This is of course 
generically impossible. 

However, it is possible to diagonalize two {K and H) of 
the three matrices simultaneously by the trick of prewhiten- 
ing, and to cross one's fingers on the third matrix (J) . The 
term prewhitening refers to the operation of multiplying a 
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signal by a function in such a way that the noise becomes 
white, or constant (Blackman & Tukey 1959 §11). Prewhi- 
tening is commonly used in the construction of Karhunen- 
Loeve modes (signal-to-noise eigenmodes), in order to al- 
low a signal and its noise to be diagonalized simultaneously 
(Vogeley & Szalay 1996; Tegmark, Taylor & Heavens 1997; 
Tegmark et al. 1998). 

Define the prewhitened reduced covariance 5Bq/3 to be 



-1/2 



(63) 



and similarly define the prewhitened 4-point and 3-point 
matrices Ma^ and La^ to be 



M = H' 



-1/2 



KH 



-1/2 



H 



-1/2 



JH 



-1/2 



(64) 
(65) 



By construction, the prewhitened 2-point matrix is the unit 
matrix, H~^^^H H~^^^ = 1. In terms of the prewhitened 4- 
point and 3-point matrices Ma/3 and La/j, the prewhitened 
reduced covariance ^ai3{ni,nj) is (compare eq. p^) 

^af3{ni,nj) = 2 [Mai3 + {ni^+nJ^)Lai3 + nY^nJ^lafi] .(66) 

The properties of the prewhitened 4-point and 3-point 
matrices M and L are examined in SM. 



3.4 FFTLog 

Several of the manipulations described in this paper involve 
transforming between real and Fourier space. Ideally, one 
would like to be able to cover several orders of magnitude in 
separation or wavenumber. The SDSS, for example, should 
be able to probe scales from 10^^ ft~^Mpc to 10^ /i~^Mpc, 
a range of 10^. If the Fourier transforms were done using 
standard Fast Fourier Transform (FFT) techniques, which 
require lineary spaced points, covering such a range would 
require 10^ points. The trouble with this is that one would 
then have to manipulate 10^ x 10^ matrices. Clearly this is 
a problem of the shoe not fitting the foot; that is, a linear 
spacing of points is not well suited to the case at hand: 
while the difference between separations of 0.01 /I'^Mpc 
and 0.02 /i~^Mpc may be significant, the difference between 
1000.01 /i"^Mpc and 1000.02 /i"^Mpc is practically irrele- 
vant. 

The problem may be solved by using an FFT method 
originally proposed by Talman (1978), that works for loga- 
rithmically spaced points, and which I have implemented 
in a code FFTLog. FFTLog is analogous to the normal 
FFT in that it gives the exact Fourier transform of a dis- 
crete sequence that is uniformly spaced and periodic in log- 
arithmic space. More generally, FFTLog yields Fast Hankel 
(= Fourier-Bessel) Transforms of arbitrary order, including 
both integral and 1/2- integral orders. FFTLog, like the nor- 
mal FFT, suffers from the usual problems of ringing (re- 
sponse to sudden steps) and aliasing (periodic folding of 
frequencies), but under appropriate circumstances and with 
suitable precautions, discussed in Appendix B, it yields re- 
liable Fourier transforms covering ranges of many orders of 
magnitude with modest numbers of points. 

Appendix B gives further details of FFTLog. The code 
may be downloaded from http://casa.colorado.edu/~ajsh 
/FFTLog/ . 



4 PREWHITENED 4-POINT AND 3-POINT 
COVARIANCE MATRICES 

4.1 Computation 

Before showing pictures, it is helpful to comment on the 
numerical computation of the 4-point and 3-point covariance 
matrices Ka/s and J^p and their prewhitened counterparts 

Mai3 and Lq/j. 

Equations ( |55[ ) and (B6|) give expressions for the 4-point 
and 3-point matrices K{ka,ki)) and J{ka,kp) in Fourier 
space, for the hierarchical model with constant hierarchi- 



cal amplitudes. These are discretized as described in §2.3 



An issue here is the calculation of the subsidiary matrix 
A{ka,ki3). This matrix Aafs is diagonal in real space with 
diagonal entries ^(ra), equation (pSl), so one way to calcu- 
late A{ka, k/s) is to start with the diagonal matrix A{ra,r0) 
in real space, and then Fourier transform it into Fourier 
space. Unfortunately the resulting Fourier transformed ma- 
trix A{ka,ki3) shows evident signs of ringing and aliasing, 
which is true whether the wavenumbers ka are linearly 
spaced (FFT) or logarithmically spaced (FFTLog). Part of 
the difficulty is that the diagonal matrix A(rQ,, r^) is liable to 
vary by several orders of magnitude along the diagonal; since 
the FFT (or FFTLog) assumes that the matrix is periodic, 
the matrix appears to have a sharp step at its boundary. 
These problems can be reduced by padding the matrix, and 
in the case of FFTLog by biasing the matrix with a suitable 
power law (see Appendix B). Still, artefacts from the FFT 
remain a concern. 

A more robust procedure, the one used in this paper, 
is to avoid FFTs altogether, and to calculate the matrix 
A{ka, kff) directly from its Fourier expression ([59[). 

A similar issue arises when prewhitening the 4-point 
and 3-point matrices K and J. The prewhitening matrix 
is again diagonal in real space, with 
diagonal entries [1 -|- ^{r)]~^^^. Thus one way to prewhi- 
ten K (say) is to start with K{ka,kf3) in Fourier space, 
Fourier transform it into real space, prewhiten M{ra, rp) = 
[1 + i{rc)]-^/''K{r^,rii)[l + e(r/3)l"'/^ and then Fourier 
transform back into Fourier space. Once again the result- 
ing matrix M{ka, kp) shows signs of ringing and aliasing. 

Again, a more robust procedure, the one used in this 
paper, is to avoid FFTs, and to calculate the prewhitening 
matrix H~^^^ = {1+A)~^^'^ directly in Fourier space. Specif- 
ically, take the Fourier expression ( ^9| ) for A{ka,ki3), add 
the unit matrix 1 to form H , and evaluate the inverse posi- 
tive square root H~^^^ via an intermediate diagonalization. 
This yields the prewhitening matrix H~^^^ in Fourier space, 
which can be used directly to prewhiten the 4-point and 3- 
point covariances matrices K ox J in Fourier space. This 
manner of constructing H^^^^ guarantees that the prewhi- 
tened 2-point covariance matrix H~^^^H H~^^^ is numeri- 
cally equal to the unit matrix 1, as it should be. Although 
this procedure is slower than using FFTs, it yields results 
that are robust with respect to range, resolution, and linear 
or logarithmic binning, and consistent with the results from 
FFTs if due care is taken with the latter. 



4.2 Prewhitened 4-point covariance matrix 

Figure 1^ shows the correlation coefficient Kap / {KaaKpfsY^^ 
(no implicit summation) of the 4-point contribution Kafs to 
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7= 1.1 r= 1.8 7= 2.9 ACDM 




Wavenumber k (hMpc ') 



Figure 4. Correlation coefficients (top) K(ka, kp)/[K{ka, ka)K{k0, kp)]^^^ of the covariance, and (bottom) M{ka, k0)/[M(ka, ka) 
M{kfj, fc^)]^/^ of the prewhitened covariance, of the power spectrum in four different models of the power spectrum. Each curve is the 
correlation coefficient at fixed kp = Ih Mpc~^, plotted as a function of ka- The three sets of panels starting from the left are for power 
law power spectra with correlation functions ^(r) = (r/5 /i~^Mpc)~''' with indices 7 = 1.1, 1.8, and 2.9, while the rightmost panel is 
for the ACDM power spectrum of Eisenstein & Hu (1998) with = 0.7, Om = 0.3, Q.f,h? = 0.02, and h = 0.65, nonlinearly evolved 
by the procedure of Peacock &; Dodds (1996). The two lines on each graph are for 4-point hierarchical amplitudes (solid) Ri, = —Ra, 
and (long-dash) R^ = Ra. Lines are dotted where the correlation coefficient is negative. The Schwarz inequality, which requires that 
the correlation coefficient be < 1, is violated by the hierarchical model with ij;, = Ra at values k <^ k' and k 2> k' . The resolution is 
128 points per decade, the same as in Figures u and pi 




Figure 2. Correlation coefficient K{ka,kfj)/yK{ka:,ka)K{kf}, 
A;^)]^/^ of the 4-point contribution K{kci,kf}) to the covari- 
ance of the power (i.e. the covariance without shot-noise) in the 
case of a power law power spectrum with correlation function 
^(r) = (r/5 /i~^Mpc)~^'*. Each line is the correlation coefficient 
for a fixed fc^, and each line peaks at ka = fc^, whereat the value 
is unity. The hierarchical amplitudes are Ra = —Rh = 4.195. The 
resolution is 128 points per decade, A log A; = 1/128. 



the (unprewhitened) reduced covariance matrix £a/3, equa- 
tion (^), for the case of a power law power spectrum having 
correlation function ^(r) = (r/5 /i~^Mpc)~^ *. Physically, 
the quantity plotted is the (correlation coefficient of) the 
covariance of estimates of power in the case of a perfect 
survey with no shot-noise, n ^ 00. 

The correlation coefflcient offers a good way to visualize 
the covariance, since a value of ( — )1 means two quantities 



Figure 3. Correlation coefficient M(ka,kfj)/[M{ka,ka)M(kfj, 
kfj)Y'^'^ of the 4-point contribution M(kc k^) to the prewhitened 
covariance of a power law power spectrum with correlation func- 
tion ^(r) = (r/5 h^^Mpc)"-*^'*. Lines are dotted where the corre- 
lation coefficient is negative. This is the same as Figure ^ except 
that the covariance is prewhitened. 

are perfectly (anti-)correlated, and the Schwarz inequality 
requires that the absolute value of the correlation coefficient 
always be less than or equal to unity. 

The Gaussian spikes evident in the curves on the left- 
ward, linear, side of Figure |^ reflect the fact that the co- 
variance of power becomes diagonal in the linear, Gaus- 
sian regime. In the nonlinear regime, the hierarchical con- 
tribution to the covariance dominates, and the covariance 
of power becomes quite broad, a point previously made by 
Meiksin & White (1999) and Scoccimarro et al. (1999). 
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It should be borne in mind that the shape of the correla- 
tion coefficient shown in Figure ^ depends on the resolution 
in wavenumber k, a point emphasized by Scoccimarro et al. 
(1999). In Figure ^ the points are logarithmically spaced 
with 128 points per decade, so Alogfc — 1/128. However, 
the correlation coefficient varies in an unsurprising way: as 
the resolution increases, the Gaussian spikes gets spikier, 
tending in principle to a Dirac delta-function in the limit of 
infinite resolution. 

Figure H shows the correlation coefficient Map/ 
(MaaMfjfj)^^ of the 4-point contribution MaS to the pre- 
whitened reduced covariance 23^/3, equations ( |63| ) and (^^, 
again for the case of a power law power spectrum having cor- 
relation function ^(r) = (r/5 /i~^Mpc)~^ *. The only differ- 
ence between this Figure and Figure H is that the covariance 
is now prewhitened. 

The prewhitened covariance M plotted in Figure ^ ap- 
pears to be remarkably narrow, certainly substantially nar- 
rower than the covariance shown in Figure ^. The Gaus- 
sian spikes again show up in the linear regime, and again 
the hierarchical contribution to the prewhitened covariance 
dominates in the nonlinear regime. The hierarchical con- 
tribution appears empirically to have a constant width of 
Ak « n/ro ~ 1 /iMpc~^, where ro — 5 ^~^Mpc is the corre- 
lation length. Thus the prewhitened covariance appears to 
become relatively narrower at large wavenumber k. 

Figure ^ shows the correlation coefficients of the covari- 
ance of the power, both straight K and prewhitened M, for 
several other power spectra. In each case the covariance of 
power with the power at A: = 1 hMpc~^ is plotted, which is 
essentially the 'worst case', where the prewhitened covari- 
ance M is relatively broadest. 

The solid lines in Figure ^ are for 4-point hierarchical 
amplitudes Rb — —Ra, while the dashed lines are for Rb — 



Ra- As discussed in §3.2, the hierarchical model violates the 
Schwarz inequality at ka 2> kp (or ka <C kp) unless Rb — 
—Ra- 

Figure M illustrates that the pattern encountered in Fig- 
ures 1^ and l^is remarkably robust over different power spec- 
tra. That is, while the covariance of the power is itself broad, 
in all cases the covariance of the prewhitened power is sub- 
stantially narrower, at least for Rb = —Ra (solid lines). Note 
that the power law power spectra illustrated in Figure ^ 
cover essentially the full range of indices, 1 < 7 < 3, allowed 
by the hierarchical model. 

The situation for Rb = Ra is muddier. Although the 
core of the prewhitened covariance is for the most part rea- 
sonably narrow also in this case, the off-diagonal covariances 
at ka ^ kfs (or ka <S fe^) are starting to become worrying 
large in several cases. Some of this behaviour is undoubtedly 
inherited from the unphysical (Schwarz-inequality- violating) 
behaviour of the ordinary covariance, and is surely not real- 
istic. Here I leave the problem with the comment that fur- 
ther investigation is clearly required, along the lines being 
pioneered by Scoccimarro et al. (1999). 



4.3 Prewhitened 3-point covariance matrix 



As discussed in §3.3, it would be ideal if the prewhitened 
3-point contribution Lap to the covariance of power were 
diagonal in the same representation as the 4-point contribu- 
tion Mal3- 
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Figure 5. Correlation coefficient Lap/ (LaaLpp)^^^ of the pre- 
whitened 3-point covariance Lap in the representation of eigen- 
functions (pa of the prewhitened 4-point covariance Map, for 
a power law power spectrum with correlation function ^(r) = 
(r/5 /i~^Mpc)~^-*. Each line is the correlation coefficient at 
a fixed nominal wavenumber kp, plotted against the nominal 
wavenumber ka, which labels the 4-point eigenfunctions ij>a or- 
dered by eigenvalue. Each line peaks at ka = kp, whereat the 
correlation coefficient is unity. The upper panel is for 4-point hi- 
erarchical amplitudes Rb = — -Ra; the lower panel is for Rb = Ra- 
Lines are dotted where the correlation coefficient is negative. The 
resolution is A log A; = 1/128. 



Figure ^ shows the correlation coefficient of the pre- 
whitened 3-point covariance Lap in the representation of 
eigenfunctions of the prewhitened 4-point covariance Map, 
for the case of ^(r) = (r/5ft"^Mpc)-^ **. The horizontal axis 
here is a nominal wavenumber ka labelling each eigenfunc- 
tion (f)a(k) of the prewhitened 4-point covariance Map- In 
practice, the eigenfunctions are simply ordered by eigen- 
value, which in most cases (see below) yields a satisfactory 
ordering by wavenumber, in the sense that the corresponding 
eigenfunctions (j)a{k) have their largest components around 

At first sight, the correlation coefficient plotted in Fig- 
ure ^ looks astonishingly diagonal at all wavenumbers, for 
both Rb — —Ra and Rb — Ra- However, as Scoccimarro 
et al. (1999) emphasize, off-diagonal elements, though they 
may be small, are many. The resolution in Figure |5| is 128 
points per decade, and the off-diagonal elements in the case 
Rb — —Ra are down at the level of ~ 1/100, which means 
that cumulative off-diagonal covariance over a decade of 
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.2 .5 1 2 5 .2 .5 1 2 5 .2 .5 1 2 5 .2 .5 1 2 5 
Nominal wavenumber (/zMpc~^) 

Figure 6. Correlation coefficient Lcp/iLaaLiSis)^^^ of the prewhitened 3-point covariance L^/j in the representation of eigenfunctions 
</)a of the prewhitened 4-point covariance M^p, at a representative nominal wavenumber kf} = l/iMpc"^. The horizontal axis is the 
nominal wavenumber ka, which labels the 4-point eigenfunctions ipa ordered by eigenvalue. The three sets of panels starting from the 
left are for power law power spectra with correlation functions ^(r) = (r/5 /i^^Mpc)""^ with indices 7 = 1.1, 1.8, and 2.9, while the 
rightmost panel is for the ACDM power spectrum of Eisenstein & Hu (1998) with Qa = 0.7, Qm = 0.3, n^h"^ = 0.02, and h = 0.65. 
Upper panels are for 4^point hierarchical amplitudes = —Ra, lower panels for Ri, = Ra- Lines are dotted where the correlation 
coefficient is negative. The resolution is Alogfc = 1/32, four times coarser than that of Figure [sj What appears to be noise in the 
curve for 7 = 2.9 results from a degeneracy of eigenvalues that mixes the correspondence between eigenfunctions 0^ and nominal 
wavcnumbcrs kn. 



wavenumber would be comparable to the diagonal variance. 
Curiously, the off-diagonal elements are somewhat smaller 
for Rf, = Ra than for Rt — —Ra- 
in ^ and thereafter the approximation will be made, 
equation (|8^), that the 3-point matrix L is indeed diago- 
nal in the representation of 4-point eigenfunctions. If L is 
not precisely diagonal, then the 'minimum variance' pair- 
weighting that emerges from assuming diagonality will not 
be precisely minimum variance. But a linear error in the 
pair-weighting will raise the variance quadratically from its 
minimum, so the pair-weighting should be close to minimum 
variance as long as L is not too far from being diagonal. In 
as discussed in i 



any case. 



7.1 



the estimate of power remains 
unbiased whatever approximations are made. 

Figure ^ shows the correlation coefficient of the prewhi- 
tened 3-point covariance Lafs in the representation of eigen- 
functions of the prewhitened 4-point covariance Map for a 
number of different power spectra, at a representative nomi- 
nal wavenumber ka — Ih Mpc~^ . The Figure illustrates that 
this correlation coefficient remains remarkably diagonal for 
all power spectra. Again, the range of power law power spec- 
tra shown covers essentially the full range 1 < 7 < 3 allowed 
by the hierarchical model. 

In the case 7 = 2.9, the off-diagonal elements of the 
correlation coefficient shown in Figure ^ appear to bounce 
around, even though taken as a whole the correlation corre- 
lation coefficient appears more diagonal in this case than any 
other. The apparent noise is caused by a near degeneracy of 
eigenvalues. Such degeneracy is not too surprising, since in 
the limit 7^8, the prewhitened 3-point and 4-point matri- 
ces Lafi and Mafs are both expected to become proportional 



to the unit matrix. Numerically, for both 3-point and 4-point 
matrices, there is a degeneracy of eigenvalues between eigen- 
functions at small and large wavenumbers (in the sense that 
eigenfunctions with nearly the same eigenvalue may have 
their largest components at either small or large wavenum- 
bers) : the eigenvalues are larger at small and large wavenum- 
ber, and go through a minimum at intermediate wavenum- 
ber. The degeneracy causes mixing of the eigenfunctions at 
small and large wavenumber, making the correspondence be- 
tween eigenvalue and nominal wavenumber ambiguous, and 
resulting in the oscillations in the off-diagonal components 
apparent in Figure ^. 



4.4 4-point and 3-point eigenvalues 

Denote the eigenvalues of the 4-point and 3-point prewhi- 
tened covariance matrices M and L by 



Ltfia — Xaipa 



(67) 
(68) 



(no implicit summation on the right hand side) so that for 
Gaussian fluctuations the eigenvalues fi^ and Aq would be 

Figure ^shows the ratio ^a/£,{ka) of the 4-point eigen- 
values /ic to the nonlinear power spectrum £,{ka), plot- 
ted as a function of the nominal wavenumber ka, which 
labels the eigenfunctions (j>a ordered by eigenvalue, for a 
power law power spectrum with correlation function ^(r) = 



(r/5/i"^Mpc) 



The eigenvalue is comparable to the 



power spectrum at all wavenumbers, /i^ ~ ^(ka)- In the 
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Figure 7. Ratio fia/Hka) of the eigenvalue /la of the 4-point 
prewhitened covariance matrix M to the nonUnear power spec- 
trum ^(ka), as a function of the nominal wavenumber ka, which 
labels the 4-point eigenfunctions </<a ordered by eigenvalue, for 
a power law power spectrum with correlation function ^(r) = 
(r/5 /i~^Mpc)~^'*. The relation between eigenvalue and nominal 
wavenumber varies with resolution. The low resolution case has 
A log A: = 1/32; the solid line is for ij^ = —Ra, the long-dashed 
line for Ri, = Ra- The high resolution case has Alogfc = 1/128; 
here the (dashed) ijj, = —Ra and (dotted) Ri, = Ra curves lie 
practically on top of each other. 



Gaussian, small regime the eigenvalue is equal to the 
power spectrum, fia = £,{ka), as expected, while in the hier- 
archical, large ka regime the eigenvalue asymptotes to close 
to 2 Ra^^ times the power spectrum, jj.a ~ 2 Ra^^£,ika)- Sim- 
ilar behaviour is found for other power spectra (not plotted), 
and for the 3-point eigenvalue Xa, which in the hierarchical 
regime asyrnptotes to Xa ~ 2Q^{ka)- 

Figure M shows the ratio Xa/ ^la of 3-point to 4-point 
eigenvalues, as a function of the nominal wavenumber ka, 
for several power spectra. Remarkably, the ratio Xa/l-i-a of 
eigenvalues is quite close to unity at all wavenumbers and for 
all power spectra. The case 7 = 2.9 is not plotted, in part 
because of the same problem of mixing of eigenfunctions 
shown in Figure ^. In any case, for 7 = 2.9 the ratio Xa/fJ-a 
differs from unity by less than 1 percent at all wavenumbers. 
Analytically, the ratio is expected to equal one in the limit 
7^3. 

In the ACDM model, the eigenfunctions (j>oi (and tpa) 
mix where the eigenvalues fia (and Xa) are degenerate, 
which happens because the ACDM power spectrum goes 
through a maximum at A: ~ 0.017 ft Mpc^^. For the pur- 
pose of plotting the ratio XaJlia for the ACDM model in 
the bottom panel of Figure ra, this mixing was avoided by 
the device of truncating the matrices Mafs and Lap at a 
wavenumber close to the peak. Mixing causes no problems 
for the evaluation of the minimum variance estimator and 
Fisher matrix of the prewhitened power spectrum in §^ and 
§1^ (so there is no need to truncate the matrices in general) , 
but mixing does muddy the physical interpretation of the 
eigenfunctions. 
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Figure 8. Ratio Xa/fJ-a of the eigenvalues of the 3-point and 4- 
point prewhitened covariance matrices, for various power spectra. 
The two lines in each case are for 4-point hierarchical amplitudes 
(solid) ijj, = —Ra, and (long-dash) R), = Ra- The horizontal 
axis is the nominal wavenumber ka , which labels the 3-point and 
4-point eigenfunctions ipa and (pa ordered by eigenvalue. The re- 
lation between eigenvalue and nominal wavenumber varies with 
resolution. The resolution is Alogfc = 1/32, except for a high 
resolution case shown for 7 = 1.8, where A log A; = 1/128. 



Curiously, the ratios iia/i(ka) and Xa/i(ka), regarded 
as functions of the nominal wavenumber ka, vary with the 
resolution A log k of the matrix, as illustrated in Figures ^ 
and|^for the case 7 = 1.8. In the Gaussian limit of small ka, 
the ratios do not change with resolution, but in the hierarchi- 
cal limit of large ka , the ratios seems to shift (to the right on 
the Figures, as the resolution increases) in such a way that 
the ratios are functions of the product fee A log k. At interme- 
diate ka, the shift is intermediate. Now the wavenumber ka 
is only a nominal wavenumber, a labelling of the eigenfunc- 
tions ordered by eigenvalue, and it is only in the Gaussian 
regime that the eigenmodes are Fourier modes and the corre- 
spondence between nominal and true wavenumber is precise. 
Still, the shift seems surprising; for example, in the limit of 
infinite resolution Alogfc — > 0, the ratio Ha/iika) plotted 
in Figure Q would shift to the right so far that fj,a/^{ka) 
would equal 1 at all finite wavenumbers. Similarly, the ratio 
Xa/fJ-a plotted in Figure ^ would shift to the right so far that 
Xa/fJ,a would equal 1 at all finite wavenumbers. Numerically, 
to the limit that I have tested it (Alogfc = 1/1024), this is 
indeed what seems to happen: both fJ,a/^{ka) and Xa/^{ka), 
hence also their ratio Xa/fJ-a, shift to the right together as 
the resolution increases, for all power spectra. 

This does not appear to be a numerical error, because 
'observable' quantities computed via the eigenfunctions (pa 
and their eigenvalues fj.a, such as the error bars attached to 
the prewhitened power spectrum X{k) in Fourier space (§^, 
appear robust against changes in resolution. 
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Figure 9. Sample of (discretized) eigcnfunctions (/ia(fc)[47rfc'^ 
A In fe/(27r)^]^/'^ of the prewhitened 4-point covariance M^p for 
a power law power spectrum with correlation function i;{r) = 
{r /5 h~^Mpc)~^'^ , and Rf, = —Ra, at resolutions of (top) 
Alogfe = 1/32 and (bottom) Alogfc = 1/128. The eigenfunc- 
tions have nominal wavenumbers ka of 0.1, 1, and lO/iMpc"^. 
In the hierarchical regime, the eigcnfunctions grow wigglier as the 
resolution increases. 



Examination of the eigenfunctions of the 4-point and 
3-point matrices M and L reveals at least part of the rea- 
son why their eigenvalues seem to shift as the resolution 
increases. Figured shows a sampling of eigenfunctions of 
the 4-point matrix M for the case ^(r) = (r/5 /i~^Mpc)~^'*, 
at two different resolutions, Alogfc = 1/32 and 1/128. 
Whereas in the Gaussian, small ka regime the eigenfunc- 
tions go over to delta-functions in Fourier space, in the 
hierarchical, large ka regime the eigenfunctions grow ever 
wigglier as the resolution increases. What seems to hap- 
pen is that, as the resolution increases, eigenfunctions at 
neighbouring nominal wavenumbers kc, strive to remain or- 
thogonal to each other, which they accomplish by becoming 
wigglier and wigglier. To the limit that I have tested it nu- 
merically, there seems to be no end to the wiggliness. Given 
that there is no asymtotic hmit to which the eigenfunctions 
appear to tend, perhaps it is not surprising that their eigen- 
values should shift systematically too. However, it would be 
nice to have a better understanding of what is going on. 



5 PREWHITENED POWER SPECTRUM 
5.1 Definition 

Given the nice properties of the prewhitened covariance of 
power established in the previous section, it makes sense 
to define a prewhitened power spectrum Xa, and a corre- 
sponding estimator Xa thereof, with the property that the 
covariance of the prewhitened power equals the prewhitened 
covariance of power. 

Define, therefore, the prewhitened power spectrum Xa 
by, in the real space representation. 



X{r 



2C(r) 



l + [l + £.{r)Y'^ 



(69) 



The expression (||) is equivalent to X{r) = 2 [1+C(r-)]^/^ -2, 
but the former expression (^^ is numerically stabler to eval- 
uate when ^(r) is small. Similarly, define an estimator Xa 
of the prewhitened power in terms of the minimum variance 
estimator ^a, equation (^), of the power spectrum by, again 
in the real space representation. 



X(r) 



2CW 



l + [l + i{r)Y'^ 



(70) 



which by construction has the property that for small 
AX(r), as should be true in the limit of a large amount 
of data, (the following equation is essentially the derivative 
ofeq. 0), 



AX(r) 



AC(r) 



[l + e(r)]V2 



(71) 



The covariance of the estimate Xa of the prewhitened 
power spectrum is given by 



(AX,AX;3 



{H-^'^a (Ai.Ais) {H-'^% = i?-,^ 



(72) 



where the Fisher matrix E""^ of the prewhitened power 
equals the prewhitened Fisher matrix of the power, equa- 
tion (H), 



E 



all 



{H''XF {H^'yi . (73) 

In it will be found convenient to deal with another 
prewhitened estimator Ya defined by 



Ya = {H' 



(74) 



The prewhitened estimator Ya has the same covariance as 

Xa 



{AYaAYp) = {\XaAXp) = E-l 



(75) 



So why not define Ya to be the prewhitened power? The 
problem with the estimator Ya is that it depends explicitly 
on the prior power spectrum ^a ■ That is, Ya in real space is 



y(r) = 



[l + ^(r)]i/2 



(76) 



which involves an estimated quantity ^(r) in the numerator 
and the prior quantity ^(r) in the denominator. Imagine 
plotting Ya on a graph. What is this quantity supposed to 
be an estimate of? Obviously Ya is an estimate of Ya = 
(Yq) = {H~^^'^)a ^13- But if one wanted to attach error bars 
to the estimate, then to be fair one should include the full 
covariance of the quantity being estimated, including the 
covariance that arises from the denominator [1 -I- ^{r)]^^^ 
in equation (|76|), not just the covariance (AYaAYis) with 
the denominator held fixed. Indeed, if one goes through the 
usual ML cycle of permitting the data to inform the prior, 
so that the estimated ^(r) is inserted into the denominator 
of equation ([t^), then it becomes abundantly evident that 
it would be correct to include covariance arising from the 
denominator. 

To avoid confusion, it should be understood that the 
quantities Ya are of course perfectly fine for carrying out 
ML estimation of parameters. In ML estimation, 'error bars 
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Figure 10. Linear power spectrum 5l(A:), nonlinear power spectrum ^(fc), and prewhitened nonlinear power spectrum X{k) for (left) 
the Qrn = 0.3 power spectrum derived from observations by Peacock (1997), and (right) the COBE-normalized ACDM power spectrum 
from the fitting formulae of Eisenstein Sz Hu (1998), with parameters as listed on the graph. The nonlinear power spectra were computed 
from the linear power spectra according to the formula of Peacock & Dodds (1996). The ACDM power spectrum is the one used in 
Figures |, |, |, |ll[ and Q 



are attached to the model, not to the data', to quote an- 
other of the refrains from the 1997 Aspen workshop on Pre- 
cision Measurement of Large Scale Structure. Whereas in 
ML parameter estimation with Xa one might form a likeli- 
hood from the 'data' quantities AX{r) — 2^(r)/{l + [1 -|- 
iir)]^^^}-'2^(.r)/{l + [l + ^{r)Y^^}, in ML parameter esti- 
mation with Ya one would instead form a likelihood from 
the 'data' quantities AY{r) = [i{r) - + ^{r)]^^'^- 

But, for the purpose of plotting quantities on a graph, 
plainly it is the prewhitened power spectrum Xa defined by 
equation ( |70| ) that should be plotted, not Ya. 

5.2 Picture 

Figure shows prewhitened nonlinear power spectra X{k), 
along with linear and nonlinear power spectra S,L{k) and 
^(fc), for the observationally derived power spectrum of Pea- 
cock (1997) with flm = 0.3, and for a ACDM model of Eisen- 
stein & Hu (1998) with observationally concordant parame- 
ters as indicated on the graph. 

The nonlinear power spectra ^(fc) were constructed from 
the linear power spectra ^^(fc) according to the formula of 
Peacock & Dodds (1996). Amongst other things, the Pea- 
cock & Dodds formula depends on the logarithmic slope of 
the linear power spectrum. Now the Eisenstein & Hu power 
spectrum contains baryonic wiggles, causing the slope to os- 
cillate substantially, whereas what Peacock & Dodds had in 
mind was a rough average slope. For the slope of the ACDM 
model in the Peacock & Dodds formula, I therefore used the 
slope of the 'no-wiggle' power spectrum provided by Eisen- 
stein & Hu as a smooth fit through the baryonic wiggles. 
The alternative of using the wiggly slope has the additional 
demerit that it amplifies baryonic wiggles in the nonlinear 
regime, which is opposite to the suppression of baryonic wig- 



gles in the nonlinear regime observed in A''-body simulations 
by Meiksin, White & Peacock (1999). 

The prewhitened power spectra X{k) shown in Fig- 
ure |l^ were computed by transforming the nonlinear power 
spectrum ^(fc) into real space using FFTLog (see Ap- 
pendix B, Fig. ^2|), constructing the prewhitened power 
X{r) from ^(r) according to equation (^), and Fourier 
transforming back. 

The prewhitened power spectra shown in Figure |l^ ap- 
pear to be interestingly close to the linear power spectra, 
X{k) « Ci(fc), another one- eyebrow- raising property of the 
prewhitened power spectrum. But surely this is just coin- 
cidence, since for a primordial power spectrum ^(fc) oc fc" 
the prewhitened correlation in the highly nonlinear regime 
should go as X(r) « 2C(r)^/^ oc r-3("+3)/2("+5) assuming 
stable clustering (Peebles 1980, eq. [73.12]), whereas the lin- 
ear power spectrum would go as r~^"'^^\ whose power law 
exponents agree only in the limiting case n —3. Still, the 
coincidence is curious. 

Figure |l^ points up one defect of the prewhitened power 
spectrum, which is that, surprisingly enough, it does not re- 
produce the linear power spectrum at the very largest scales 
(small k). Indeed the prewhitened power goes negative in 
the Peacock (1997) case a.t k ^ 0.0023 /i Mpc"\ and in the 
ACDM case at k ^ 0.00021 /i Mpc~\ This turns out to be 
a generic feature of the prewhitened power spectrum if the 
true power spectrum goes to zero at zero wavenumber, as 
is true for Harrison-Zel'dovich models, ^(fc) oc fc as fc — > 0. 
For if it is true that the power spectrum ^(fc) goes to zero 
at zero wavenumber k 



lim5(fc)= / 5(r)47rr^dr = (77) 
Jo 
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then it follows that the prewhitened power must go to a 
negative constant at zero wavenumber 



Now {H~^^^)2 commutes with Dl^ Ofj, since both are simul- 
taneously diagonal in real space. It follows that the Fisher 
matrix E = H^^^ F H^^^ of the prewhitened power, equa- 
tion ([73|), is, in the FKP approximation, 



lim X(k) 

k-,0 



2C(r) 



l + [l + e(r)]i/2 



47rr dr < 



(78) 



since the factor 2/{l + [1 -|- ^(j')]^^'^} in the integrand is 
less than one for all positive (,{r), and greater than one for 
all negative ^(r). It is not clear what to do about this, if 
indeed anything needs to be done. Adding a constant to 
X{k) and X(k) (which would leave AX, hence the covari- 
ance (AX AX), unchanged) would spoil the nice behaviour 
of the prewhitened power in the nonlinear regime. 



6 FISHER MATRIX OF PREWHITENED 
NONLINEAR POWER IN A SURVEY 

It was found in §^ that the prewhitened reduced covariance 
*B of power appears to have some unexpectedly pleasant 
properties: first, the prewhitened covariance is surprisingly 
narrow in Fourier space; second, the 4-point and 3-point 
contributions M and L, equation (^), to the prewhitened 
reduced covariance 25 are almost simultaneously diagonal 
(the 2-point contribution is by construction the unit matrix, 
so is automatically diagonal in any representation); third, 
the 4-point and 3-point eig envalues and Aa, as defined 
by equations (^^ and (^^, are approximately equal; and 
fourth, all these results hold for all power spectra tested. 

It should be emphasized that the pleasant properties 
of the prewhitened power are not perfect, and that they 
are premised on the validity of the hierarchical model with 
constant hierarchical amplitudes, which as discussed in !;4.2 
is certainly wrong at some level. 

These pretty properties lead to an approximate expres- 
sion, equation (p3|), for the Fisher matrix of the prewhitened 
nonlinear power spectrum of a galaxy survey, which looks 
the same as the FKP approximation to the Fisher matrix of 
the power in the linear, Gaussian case, with the difference 
that the eigenmodes of the prewhitened covariance M of the 
nonlinear power take the place of the Fourier modes in the 
linear case. 



6.1 Fisher matrix 

To the extent that the prewhitened 4-point and 3-point ma- 
trices M and L are simultaneously diagonal, the prewhi- 
tened reduced covariance matrix '^^iiijii^nj) is diagonal in 
the representation of eigenfunctions (f)a of M and L, with 

'Ba,3(ni,nj) « 21q^ [/i^ -I- (n7VnJ^)Ac, -l-nr^wj^] . (79) 

To the further extent that \a ^ fJ,a, the prewhitened covari- 
ance matrix ^c,^{ni, fij) is just 



^Ba/3(ni, Tij) 2 lafs (/ia + rij ^)(^iQ -I- ^) 



(80) 



The Fisher matrix F"^ of the power spectrum is given 
in the FKP approximation by equation (^). In terms of 
the prewhitened reduced covariance Q5a^, the Fisher matrix 
F"'' is 
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(81) 



E 



(82) 



Like F""^ , the prewhitened Fisher matrix E"^ is asymmet- 
ric, inheriting its asymmetry from the FKP approximation, 
equation (^). 

To the extent that the approximation (|^) to <B is true, 
it follows from equation (fe2) that the Fisher matrix Eap 
of prewhitened power in the FKP approximation is, in the 
representation of eigenfunctions (j)a of the prewhitened co- 



variance, 



(83) 



Ea0 ^ / (l>a{r)(f)i3{r)R{r;iia)4Tvrdr 
Jo 

where R(r;^a) are FKP- weighted pair integrals (commonly 
denoted {RR) in the literature, for Random-Random) 



S:iD{rij - r) 



2[^Q +n(rO-i][^„ -f fi(rj)-i] 



d^r^dVj (84) 



the integration being taken over all pairs of volume elements 



ij separated by 



r in the survey. 



The FKP approximation to the Fisher matrix Eafi of 
prewhitened power, equation (^), takes the same form as 
the FKP approximation to the Fisher matrix of the power 
spectrum for Gaussian fluctuations derived in §5 of Paper 1 
and computed in §3 of Paper 2. The difference is that the 
eigenfunctions 4>a{r) and their eigenvalues /1q here take the 
place of the Fourier eigenfunctions jo{kar) and their eigen- 
values (,{ka) in the Gaussian case. 

6.2 Numerics 

Equation ( ^3| ) for Eci3 involves the eigenfunctions (f>ci{r) of 
t he p rewhitened 4-point matrix AI in real space, whereas in 
§4.1 it was suggested that the most robust way to compute 
M is in Fourier space. The problem is that FFTing the ma- 
trix Af from Fourier into real space is liable to introduce 
ringing and aliasing, which one would like to avoid. 

A more robust procedure is not to FFT M into real 
space, but rather to FFT the pair integrals R{r; fj,a) into 
Fourier space; this is the same procedure adopted in §3 of 
Paper 2 (except that R here is 1/2 that of Paper 2). If /1q is 
treated, temporarily, as a constant, then equation (^) can 
be transformed into real space to yield the diagonal matrix 

E{r,r';fia) = Ssoir - r') R{r; fia) . (85) 

Beware of equation (jsH])! It does not signify that the 
Fisher matrix is diagonal in real space, because the con- 
stant jj,a is different for each row of the Fisher matrix 
Eap. The Fourier transform of E{r,r'; fia) is E{k, k' ; ^a) = 
J jo{kr)jo{k'r) R{r; fia) iirr^dr, which simplifies to 



E{k, k'; fia) = ^ [R{k-k'; ^a) - R{k + k'; ^a)] 



(86) 



where R{k; /la) is the 1-dimensional cosine transform of 
R{r; fj-a) 



R{k; fia) = 2 I cos{kr) R{r; fia) dr 
Jo 



(87) 
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Transforming E{k, k' ; fia) into ^a-space gives 



Eal3 = / (t>a{k)(t)i3{k') E{k,k' ; ^a) 



A-rvk'^dkA-Rk'^dk' 

The cosine transform R(k;^a), equation can be done 

with either FFT or FFTLog; both work welL To ensure that 
R{k; fia) remains accurate at large (and small) wavenum- 
bers k, it helps to extrapolate R{r; fia) to small (and large) 
separations r before transforming. The transformation into 
(f>a space, equation (|88|), is done by discrete summations. 

Evaluating the Fisher matrix Ea/s with equations (p^- 
( pij ) successfully eliminates ringing and aliasing, but it intro- 
duces another problem. The problem is that equation ( |86| ) 
is liable to overestimate the value of E{k, k' ; fia) along the 
diagonal = fc' if the gridding in fc-space is too coarse to re- 
solve the diagonal properly, as typically occurs at moderate 
and large k with logarithmic gridding. What is important is 
that the integral of E{k, k'; fi^) over the diagonal be correct. 
Integrating E{k, k'; fia) over k' yields 



E{k, k'\ fj,a) k'dk' 



2tt 



R{k'; fia) dk' 



(89) 



The integral on the right can be done conveniently and re- 
liably by sine transforming (with FFT or FFTLog) the pair 
integral 



R{k'; Ha) dk' 



sm(fer) — ^ — - — - dr 



(90) 



Discretized ( |2.3D on a logarithmic grid of wavenumbers 
fc, the continuous matrix E{k, k' ; /la) becomes Ej.^/(ua,) = 
E{k,k';fia)4:n{kk'f^^Alnk/(2nf, and equation (||) be- 
comes 



Y^{k'/k)^^^Eky{fic) ^ ^ I R[k';yia)dk' 



(91) 



Numerically, if the left hand side of equation (pl|), with 
Efcfc'(/Xa) discretized from equation (js^), exceeds the right 
hand side of equation (^ij) , evaluated by equation ( |90| ) , then 
the value of the diagonal element Ekkipa) should be reduced 
so that the sum is satisfied. Ultimately, this procedure yields 
error bars on decorrelated band-powers (Paper 4) that are 
robust with respect to range, resolution, and linear or loga- 
rithmic binning. 



6.3 Coarse gridding 

Typically the pair integral R{r; Ha) is broad in real space, 
so its cosine transform R{k; fj.a) is a narrow window about 
fc ~ with a width comparable to the inverse scale length 
of the survey. It follows that the matrix E{k,k'; fj,a) given 
by equation (^^ is likewise narrow in fc-space, with a width 
comparable to the inverse scale length of the survey. More- 
over the sum in equation ( |9l| ) approximates R(0; Ha) at 
wavenumbers exceeding the inverse scale length of the sur- 
vey, which is to say at all except the largest accessible wave- 
lengths: 

^{k'/k)^^^Ekk'{Hc)^R{0;fJ.<^) for fc>scale-^ (92) 
k' 

where _R(0;/1q) is the pair integral at zero separation 



R{0;iic) = 



d\ 



2[Ha + fi{r)-^] 



(93) 



Thus if the matrix E^k' (Ha) is discretized on a grid that 
is coarse compared to the inverse scale length of the survey, 
then it is approximately proportional to the unit matrix 



Ekk' (Mo 



R{0;Ho 



The resulting discrete Fisher matrix Ea/s, equation 
diagonal in the (/>a,-representation 

Eaf) ~ ial3 R{0', fJ-a) ■ 



(94) 



(95) 



The result ( p5[ ) is analogous to that obtained by FKP for 
Gaussian fluctuations. 

Of course if this diagonal Fisher matrix, equation (p5[), 
is transformed back into Fourier space, then it is no longer 
diagonal. That is, equation ( p5|) asserts that the Fisher ma- 
trix of the prewhitened nonlinear power spectrum is approx- 
imately diagonal in ^^-space, not in Fourier space. 



7 ESTIMATE OF PREWHITENED 

NONLINEAR POWER IN A SURVEY 

T.l Unbiased estimate 

'In the case of a Gaussian distribution... rather than remov- 
ing the bias we should approximately double it, in order to 
minimize the mean square sampling error' - E. T. Jaynes 
(1996, sentence containing eq. 17-13). 

It is convenient to start out by considering the prewhi- 
tened estimator Ya defined by equation (|7^). The minimum 
variance estimator ^o, of the power s pec trum in the FKP 
approximation is given by equation (M). Translating this 
equation into prewhitened quantities, one concludes that the 
minimum variance prewhitened estimator Ya in the FKP 
approximation is, in terms of the prewhitened reduced co- 
variance *B, equation (p3|), and its associated Fisher matrix 
E, equation (|73|), 



(96) 



The estimator Ya is minimum variance if and only if 
is minimum variance, since Ya = {H~^^^)a ^f3 is a linear 
combination of (^a- 

Now the estimator Ya, equation (|9^), is intended to be 
an estimate of Ya = iH~^^^)a £,i3- But is that really true, 
given the various approximations? It will be true provided 
that the estimator is unbiased, meaning that the expectation 
value of the estimator is equal to the true value 



Ya 



(97) 



The expectation value of the estimator Ya given by equa- 
tion (SE) is, since {5i5j — Nij) — -D^C" according to equa- 
tion (I3), 

{Ya) = E-',^~'''-'{H~'^y,DlWi^^, 



E-l^-^^'-'DljDl.Y, 



(98) 



where the second line follows because {H ^''^)^ commutes 
with D'l^, both being diagonal in real space. It follows 
that the estimator Ya will be unbiased, {Ya) ~ Ya, provided 
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that the Fisher matrix E"'^ is taken to satisfy the asymmet- 
ric equation (p2[), not, for example, a symmetrized version 
of that equation. 

An important point to recognize here is that an esti- 
mate Ya of the form ( [96[ ) will be unbiased for any a priori 
choice of the matrix 58, regardless of the choice of prior 
power regardless of the hierarchical model, regardless 

of the FKP approximation, and regardless of the approxima- 
tion (such as eq. jsoj) to 5B, just so long as the matrix E in 
the estimator is interpreted as satisfying the unsymmetrized 
equation (p^). Ultimately this property of being unbiased is 
inherited from the basic prior assumption that galaxies con- 
stitute a random, Poisson sampling of an underlying sta- 
tistically homogeneous, isotropic density field, so that the 
product of overdensities SiSj at any pair of points ij sepa- 
rated by provides an unbiased estimate of the correlation 
function ^(r^). Note that the presumption here is that the 
galaxies sampled are an unbiased tracer of the galaxy den- 
sity itself, not necessarily of the mass density. 

Interpreting the estimator Ya , equation (bd) , as involv- 
ing the asymmetric matrix E°'^ , equation (|8^, should be 
regarded not as changing the estimator to make it unbi- 
ased, but rather as interpreting the estimator correctly. If 
instead the estimator Ya were interpreted as involving the 
symmetrized Fisher matrix E^°'^^ = Sym(-Q,^)iJ"^, for exam- 
ple, then the expectation value of the estimator would be 
{Ya) = E-J-^^E'^'''Yj, which is not the same as Ya, although 
of course it should be almost the same to the extent that 
Eai3 is almost symmetric. 

It is convenient to introduce yet another estimator Z" 
related to the estimator Ya by 



7.2 Numerics 



Ya = E-'^Z^ 



As in yi.2| , to avoid potential problems of ringing and alias- 
ing, i t is p robably better to evaluate the estimator Za, equa- 
tion (101), by means of an expression that involves the eigen- 



functions (j>a{k) in Fourier space rather than the eigenfunc- 
tions (jyaij) in real space. 

If ^a is treate d, te mporarily, as a constant, then trans- 
forming equation (101) into real space yields 



The Fourier transform of this is 

n ' 2 

Z{k-iia) = / jo{kr) Z{r- /la) i-^r dr 

Jo 



(103) 



(104) 



in terms of which the estimator Za, equation (101), is 
47rfc^dfc 



Za — 



0a(fc) Z{k; ^a) ■ 



(2^) 



(105) 



space, equation (lOE), is done 



The transformation into 

by discrete summation. 

The advantage of e quation (105) over t he n ominally 
equivalent equation (101) is that in equation (105) it is the 
data that are Fourier transformed, Z{r ;^ia) —> Z{k;fia), 



(99) 



equation ( 104), whereas in equation (101 ) it is the eigenfunc- 
tions of the matrix M that must be transformed, 4'a{k) — > 
(f)a{r). While the two methods would yield identical results 
for Za if the same unitary Fourier transform were applied 
in both cases, in reality it may be advantageous to have 
the freedom to Fourier transform the data the best way one 
can, without regard to the irrelevant question of how the 
eigenfunctions (j)a behave when Fourier transformed. 



In the FKP approximation, the estimator Z" is 

= ^-'^°'''{H''^'ypDl^{5^5J - N,,) . (100) 

If the approximation ( pO|) t o the prewhitened covariance 
*B is used in the estimate nTOOl) of Z, then, in the represen- 
tation of eigenfunctions ^a. 



Za 



(f)a{r)S{r; fia 
[l + ^(r)]i/= 



■ inr dr 



(101) 



where S{r;fia) is the FKP-weighted integral over pairs of 
overdensities S{ri)5{rj) at points ij separated by rij = \ri — 
rj\ = r (commonly denoted {DD) - 2{DR) + {RR} in the 
literature, D for data, R for random) 

S{r;^a) ■- 



Sso{r,-r)Sir.)Sir) , 
[fia + n{ri)-^[Ha+fi{rj)-^ ' 



The shot-noise Nij is excluded from equation (102) by ex- 
cluding from the integration the contribution from self-pairs 
of galaxies, which of course have zero separation. The as- 
sociated asymmetric Fisher matrix E"'^ is given by equa- 
tion (H). 

Equation (1102) is expressed as an integral over pairs 
of overdensities S(ri)5{rj) in real space. One could just 
as well express S as an integral over pairs of overdensi- 
ties 5{ki)5{kj) in Fourier space, or pairs of overdensities 
5{ki,£i,mi)S{kj,£j,mj) in spherical harmonic space, if one 
found it more convenient. 



7.3 The covariance of Z" 

It will now be argued that the covariances of the estima- 
tors Z" and Y°' axe approximately equal to, respectively, 
the sy rnme trized Fishe r matrix E'^°'^\ and its inverse, equa- 
tions ( 109 ) and (115). It seems worthwhile to go through 
the arguments rather carefully. As a general rule, one should 
estimate error bars as accurately as possible; but if some ap- 
proximation is necessary, then one would prefer to err on the 
conservative si de o f overestimating the true errors. 

Equation (109) will now be derived, commentary on the 
derivation being deferred to th e end . The covariance of the 
estimate Z" is, from equation (100), 



(106) 



\Loi\), while e:;^.fc, 




in which 53 is the approximate prewhitened reduced covari- 
ance matrix ( |80[ ) used to construct the estimate Z" , equa- 
tion 

(g 
("l| 



^'^^ is the true covariance matrix, equa- 
To the extent that the FKP approximation, equa- 



ls valid for equation (106) reduces to 



{AZ°'AZ 



-1/2NC 



1^7/ — — 



CqFKP/- 



FKP/ - - 
7€ 



05 



(107) 
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where is the FKP covariance, equation (^^, and 



H ^^^([^^^ H ^^"^ is its prewhitened counterpart, 



equation (pq). Note that going from equation ( [LOq ) to the 
second expression in equation (107) includes, as part of 
the FKP approximation, the assumption that fih and fii in 
S~^'''^(nfc, fij) are approximately co nsta nt. The expressions 
on the right hand side of equation (107) are not symmetric 



in a/9, because of the asymmetry of the FKP approxima- 
tion (^). To the further extent that the prewhitened co- 
variance *B^^^ equals the approximation 03, equation (po|), 
the covariance (A^" AZ^Weduces to the asymmetric ma- 
trix E given by equation (p3) 

(AZ^AZ") = E""!^ (108) 



the asymmetry of the right hand side being inherited from 
the FKP approximation. An equally good approximation 
to the covariance would be the same expression (108) with 
the indices swapped on the right hand side, a /3. 
Thus it seems reasonable to conclude that the covariance 
(AZ^AZ'^) should be approximately equal to the the sym- 
metrized Fisher matrix E^°'^^ 



Sym E 

(all) 



(109) 



Several comments can be made about the accuracy of 
the approximations made in the above derivation. 

Firstly, one partial test of the validity of the FKP ap- 
proximation is the degree of asymmetry of the asymmetric 
Fisher matrix E°'^ , equation (p3). If the survey is broad in 
real space, which is the condition for the FKP approximation 
to hold, then the pair integral R(r\ Ha) in the integrand on 
the right hand side of equation (B3) will be a slowly varying 
function of pair separation r, so that the matrix E"'^ will be 
nearly diagonal, hence symmetric. The test is not definitive 
because E°'^ would be symmetric in any case if fia = But 
in practice /Xa ~ ?(fca) in both linear and nonlinear regimes. 
Figure ^ and realistically the power spectrum ^(fc) varies 
substantially, so the consistency test should be indicative. 

Secondly, one of the weaknesses of the FKP approxi- 
mation is that it fails to deal with sharp edges - as typ- 
ically occur at the angular boundaries of a survey - cor- 
rectly. The FKP approximation tends to overestimate the 
variance contributed by regions near boundaries, since it 
assumes that those regions are accompanied by more corre- 
lated neighbours than is actually the case. Thus, at least as 
regar ds ed ge effects, the FKP approximate covariance, equa- 
tion (107), should tend to overestimate the exact covariance. 



equation (106), of the approximate estimate Z". 

Thirdly, it is possible to check the accur acy of the ap- 
proximation made in going from equation (107) to equa- 
tion (108). The approximation involves setting *B ^^05^^ = 
1, whereas comparing equation (||) for S^^^ to the approx- 
imation ( po| ) for 05 shows that this quantity is in fact, in the 
representation of eigenfunctions of the 4-point matrix M, 



{vvFKP/- 



1" + 



{pLa + )(fla + ) 



(110) 



The correction term on the right hand side of equation ( 110 ) 
should be small to the extent that the 3-point matrix Lap is 
near diagonal in this 4-point representatioii, with eigenvalues 
\a ~ fJ,a, as was found to be the case in 



If desired, one coul d u se the expression on the right 
hand side of equation (IIC) to compute a more accurate 

based on equa- 
However, if one 



approximation to the covariance of Z°' , 
tion (107) rather than on equation (10^) 
were willing to go to the trouble of computing a correction 
from equation (IIC), then one would probably be willing 
to revert to equation (IOC), and to integrate 03~^"''(ni, fij) 



numerically over all pairs of volume elements ij in the sur- 
vey, inverting *B numerically for each pair fii, nj of values 
of the selection function. This latter procedure is in fact the 



gourmet recipe of §9.1 



Fourthly, the approximation \a ~ fJ.a adopted in the 
approximation ( ^o| ) to 05 tends to overestimate the true 
eigenvalues of the 3-point matrix L, according to Fig- 
ure ^ This should lead to a slight overestimate of the vari- 
ance. In the realistic ACDM case, Figure ^, the approxima- 
tion Xa fj,a overestimates the true eigenvalues by at 
worst 20 percent, at moderately nonlinear wavenumbers k. 
This 20 percent overestimate is diluted to at worst 10 per- 
cent because the 3-point variance contributes at most half of 
the combined 2-point, 3-point, and 4-point variance, where 
the selection function satisfies = fj^a- The overestimate 
is further diluted because in practice the selection function 
varies, and is unlikely to sit everywhere near the worst value. 

The conclusion is that t he c ovar ianc e of the approxi- 
mate estimator Za, equati on ( 101) or (105), should be given 
approximately, equation (|l09|) , by the symmetrized Fisher 
matrix E^°'^^ of equation (|83|), and that if anything this co- 
variance is likely to be on the conservative side of the true 



covariance. 



7.4 The covariance of 



From the expression ( |10£| ) for the covariance of Z", one 
might conclude (falsely) that the covariance of Y^, equa- 
tion (p9[), is 

{AY^AYp} = E-.] E^-^'~> E-pl . (Ill) 

A more direct derivation of the covariance of Ya, along the 
lines of equations (106)-(109), lead s to the same (false) con- 
clusion. The analogue of equation (108) is 



{AY^AYp} = E^^ 



(112) 



with the asymmetric matrix E on the right hand side. At 
this point one might be inclined to symmetrize this equa- 
tion ( pT^ ), as was done for {AZ^AZ^^) in equation ( p^ , 
writing 



(Ay,A-f»= SymE-^ 



(a/3) 



(113) 



The symmetrized inverse Sym^^^jiS^^ of the asymmetric 
Fisher matrix is to be distinguished from the inverse E'^afs) 
of the symmetrized Fisher matrix. But it is not hard to show 
that 



Sym E-l = E^}^ E^'''^ E^l 

(aP) 



(114) 



Thus equations (111) and (Hi) are identical. However, both 
equations are wrong. 

The problem is that, while the Fisher matrix E remains 
well-behaved in the presence of loud noise, with near zero 
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eigenvalues, its inverse becomes almost singular. Con- 
sider the example of some noisy mode, for which the eigen- 
value of the Fisher matrix is almost zero. It may well happen 
that the asymmetric Fisher matrix E"''^ is numerically non- 
singular, but that, because of approximations or numerics, 
the computed eigenvalue of the s ymm etrized Fisher matrix 
S'"''^ is exactly zero. Equation (lllll) would then say that 



the variance of the noisy mode is zero (for if the determinant 
of the symmetrized Fisher matrix is zero, \E^°''^'> \ = 0, while 
the determinant of the asymmetric Fisher matrix is fi nite 



|^a/3| then the determinant of the variance in eq. [ ill 
is zero). This is plainly absurd. 

It is safer to take the covariance of Ya to be approxi- 
mately equal to the inverse of the symmetrized Fisher matrix 



(Ay^Af^) = E-\ 



(115) 



Here a noisy mode will always reveal itself by its small eigen- 
value. 



is obscure, and in the second place because the eigenfunc- 
tions can mix where their eigenvalues are degenerate. 
Since ~ £,{ka), such mixing in practice occurs between 
wavenumbers ka where the power ^(fc^) is the same, which 
happens to either side of the peak in the power spectrum. 
Figure ^ Perhaps in the future a better understanding of 
the eigenfunctions (j)a will emerge, amongst other things al- 
lowing mixing to avoided, but in the meantime these prob- 
lems remain. 

The natural solution is to decorrelate the prewhitened 
power X{k) in Fourier space. As seen in §^ the covariance 
of the prewhitened power is encouragingly narrow in Fourier 
space, narrow enough that the decorrelation matrices will be 
narrow, so that the decorrelated band-powers can be inter- 
preted as estimates of the prewhitened power over narrow 
intervals of wavenumber k. In contrast to the prewhitened 
power Xa in the ^^-representation, the prewhitened power 
X{k) in Fourier space has a clear interpretation, and there 
is no problem arising from mixing of eigenfunctions. 



7.5 Convert to Xa 

For the purpose of constructing uncorrelated quantities to 
be plotted on a graph, it is desirable to compute the pre- 
whitened power spectrum Xa- 

To co mpu te Xa, start from the estimate Za given by 
equation (105), transform this into Ya = E~p , equa- 
tion (B9|), thence into the power spectrum C« = {H-^'Yc^y^, 
equation (^) , and thence into the prewhitened power spec- 
trum Xa, equation ([to]). 

The covariance of the prewhitened powe r Xa is, by con- 
struction, the same as that of Ya, equation (115), 



{AXa^Xp) = 



/3) 



(116) 



the inverse of the symmetrized Fisher matrix of the prewhi- 
tened power. 

The estimator Xa of prewhitened power, equation ([to]), 
is a nonlinear transformation of the estimator of power, 
and is therefore biased if is unbiased. However, the es- 
timator Xa is unbiased in the asymptotic limit of a large 
quantity of data. 



T.6 Decorrelate 

One final step remains, which is to process the measured 
prewhitened power spectrum Xa into a set of decorrelated 
band-powers. How to accomplish such decorrelation is de- 
scribed in Paper 4. 

One possibility would be to decorrelate the power spec- 
trum ^(fc) itself. This is a bad idea, because the power spec- 
trum is highly correlated in the nonlinear regime, so the 
decorrelation matrices would be broad, with large negative 
off-diagonal entries, making it impossible to interpret the 
decorrelated band-powers as representing the power spec- 
trum over some narrow band. 

Another possibility would be to decorrelate the prewhi- 
tened power Xa not in Fourier space but rather in the rep- 
resentation of eigenfunctions (f)a of the prewhitened 4-point 
matrix M. Again this seems not so good an idea, in the first 
place because the physical meaning of this representation 



8 THE FULL FKP 

Sections ^ and ^ invoked not only the FKP approximation, 
but also the simplifying approximation ( |80| ) to the prewhi- 
tened reduced covariance 05. How much more work would it 
take to compute the minimum variance estimator and Fisher 
matrix of nonlinear power making only the FKP approxima- 
tion and no other approximation? The question is of both 
didactic and practical interest. 



8.1 Fisher matrix 

The FKP approximation to the Fisher matrix of the power 
spectrum, equation (psj), looks simplest expressed in real 
space: 

F{ra,ri3) = j ^^^{ra,ri3-,ni,nj)5-iD{rij - rp) A\iA\j -(117) 

The corresponding expression for the FKP approximation 
to the Fisher matrix of the prewhitened power spectrum, 
equation (^), is 

r 

E{ra,rf3) = / ^^^{ra,ri3;nt,nj)S3D{rij-rfj) d\id\j .(118) 



These are 5-dimensional (thanks to the Dirac delta-function) 
integrals over pairs of volume elements ij separated by 
Tij = \ri — Vjl — r/3 in the survey. The integrals are ac- 
tually quite doable. If, as is typical, the selection function 
n(r) separates into the product of an angular mask and a 
radial selection function, then the 3-dimensional angular in- 
tegrals can be done analytically (Hamilton 1993), leaving 
a double integral of €~\ra, r^; fii, fij) or '^~\ra, rp; fii, Uj) 
over the radial directions. The matrices (L{ra,r0;ni,nj) or 
^ {r a, rf3; fii, n j ), disaetized (§2.3) over a grid of separations 
Va and rfj, can be inverted numerically for each pair of values 
of the selection functions fii and n, . 



The problem with equations (117) or (118) is that expe- 
rience (§ g4.l| , 3.2) suggests that discretization of the matrix 
^{ra,ri3;fii,n,j) or Q5(rQ,, r^g; n^, %) in real space is liable to 
introduce ringing and aliasing in Fourier space, defeating the 
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aim of constructing an accurate Fisher matrix of the power 
spectrum. 

A possibly more robust procedure would be to follow 
more closely the program described in §^ and §^ In Fourier 
space, the FKP approximation to the prewhitened Fisher 
matrix, equation (B2), is 



E{ka,kfi) = 



*B \ka, k-,; Hi , nj)jo{k~,rij)jo{ki3rij) d^r^d^' 



(119) 



ink^^dk-y 
' (27r)3 



This integral might be evaluated as follows (since I have not 
actually carried through this program, I cannot say for sure 
that it would work without a hitch). Firstly, compute the 
matrix of pair integrals 



R{r; ka,kp 



= J ^-\k^,kis; 



Hi , nj)S3D (nj-r) d^r^dVj (120) 



for many pair separati ons rg — r. These pair integrals 
R{r;ka,ki3), equation (120), are analogous to the FKP- 
weighted pair integrals R{r; ^a), equation (^^. Next, cosine 
transform (e.g. with FFTLog) the pair integrals 



R{k- ka,kp) 



cos{kr)R{r; ka, kp) dr 



(121) 



analogously to R{k; ^c,), equation (p7|). Finally, compute the 
prewhitened Fisher matrix Eikckp) by integrating 

E{k^,k0)^ (122) 
1 f°° ~ 

) - R{k^ + ki3; k 

In practice, the matrix 55(fca, fe/j; fii, nj) in equa- 
tion ( |l2C| ) must be inverted on a discrete grid of wavenum- 
bers k. Similarly, the integral over fc-y in equation (122) 



should be done a s a d iscrete sum. Specifically, if the matrices 
are discretized (§2.3) on a logarithmic grid of wavenumbers, 
so that R{k; ka,k^ discretizes to Rk^kfiik) — R{k; ka, kjs) 
47r(fcc,fc/3)^/^Alnfc/(27r)^, and the Fisher matrix E{ka,kp) 
discretizes to Ek„k n = E{ka, kfi) An^kakp)^^^ Alnk/{2Tvy\ 
then equation (|l22|) becomes 



^^^y^ [Rk^k-,{k.-kf,) - Rk^k,{k., + kp)] .(123) 



It may be anticipated that, as in §3.2, equation (123) 
will tend to overestimate the diagonal elements E^^fe^ if 
the gridding of the matrix is too coarse to resolve the di- 
agonal properly. Int egrat ing the continuous Fisher matrix 
E{ka,k0), equation (1122), over fc^ yields 

R{k;ka,kj)dkkjdk^ .{124) 



1 



E{ka,kp)ki3dki3 = - 

'^"'0 "'0 



Discretized, equation (|12J) becomes 



TV ^ \kj 



1 . 1/2 fk-i 

' Rfc„fe,(fc)dfc .(125) 





The integral over k on the right hand side of equation (125) 
can be done conveniently as a sine transform (e.g. with 
FFTLog) of the pair integral 

J Rk^k-,(k)dk = 2 j sin(fcr)-^^i^dr . (126) 



If the sum on the left hand side of equation (125) exceeds 
the right hand side, then reduce the diagonal element E^^fe^ 
so that the sum is satisfied. It is fine to evaluate the sum on 



the right hand side of equation (125) as a discrete sum over 
k^, rather than as a continuous integral, because Rka,kjyk), 
equation (120), inherits its behaviour from '^k^k^y, which, 

from the 4-point and 3-point 



if constructed, equation 

matrices M and L as discussed in §^^, should behave cor- 
rectly near the diagonal even if the resolution is too coarse 
to resolve the diagonal. 



8.2 Estimate of power 

The FKP approximation to the minimum variance estimator 
of the power spectrum, equation (^^, again looks simplest 
when expressed in real space: 



with 

T{ra) = / €^\ra,rij\ni,nj)5{ri)5{rj)d\id\j 



(127) 



(128) 



As usual, the prewhitened estimator Y = H~^^^^, equa- 
tion (174), is related to the estimator Z hy Y = E~^Z, equa- 
tion (M). The FKP approximation to the estimator Z is, 
equation (100), 



Equations (128) and (129) are 6-dimensional integrals over 
pairs of volume elements ij in the survey. But once again 
one may anticipate that discretization of the matrices 
^{ra,rp;ni,nj) or 58(rc«, r^j; n;, fij) in real space would in- 
troduce ringing and aliasing in Fourier space, defeating the 
aim of constructing an accurate estimator of the power spec- 
trum. 

Again, it seems likely that it would be more robust 
to work with prewhitened quantities in Fourier space. In 
Fourier spac e, th e FKP approximation to the estimator Za 
is, equation (100), 



Z{ka) 



(130) 



05 \ka, kp; n„nj)jo{k0U,)5{r ,)S{rj) ^3 3 47rfc^dA:^ 

Q TiCi Vj " 



(2^) 



One way to evaluate this integral might be as follows. Firstly, 
compute the matrix S{r;ka,kis) of integrals over pairs of 
overdensities 5{ri)5{rj) at many separations rij — r 



S{r; kck/i) = 



(131) 



nj)S3D{rij-'r)S{r^)S{rj) dV^d^j 



which may be compared to equation (102). Next, prewhiten 
(compare eq. |10J]) 



Z{r; ka,kfj) 



S{r; kg, kp) 
[l + C(r)]i/2 



(132) 



and Fo urier transform, e.g. with FFTLog, (compare 
eq. jioi) 



Z{k-ka,kp)^ I jo{kr)Z{r;ka,kp)-iTvr^dr . (133) 
/o 
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Actually it suffices to do this Fourier transfo rm f or k — kp 
only. Finally, the estimator Z{ka), equation (130), is 



Z{k^) = / Z{kp-k^,kp) 



AiikjjAkfj 
(27r)3 



(134) 



The integral in equation (134) should be done as a discrete 
sum. If discretized (§2.3) on a logarithmic grid of wavenum- 
bers, so that Z{kct) discretizes to Zfc^ = Z(feQ)[47rfc„ 
Alnfc/(27r)^]^/^ and Z(k;ka,kj3) discretizes to Zk^kJk) = 



Z(fc; fee, fc^)47r(fcc,fc^)^''^Alnfc/(27r)^, then equation ( |l34| ) is 

1/2 



inkpAlnk 
(27r)3 



(135) 



9 RECIPES 

This section summarizes the results of previous sections 
into logical sequences of practical steps needed to estimate 
the prewhitened nonlinear power spectrum from an actual 
galaxy survey. The end result is a set of uncorrelated pre- 
whitened nonlinear band-powers with error bars, over some 
prescribed grid of wavenumbers k. 



There are three versions of the recipe, gourmet (§9.1), 
fine (Sp^, and fastfood (§U). AU the methods use the FKP 



approximation, equation {iw. Thus one should imagine that 
there is also a haute-cuisine method, which might be brute- 
force, or it might be some clever procedure that apodizes 
edges. 

First a disclaimer. The methods described herein do not 
take into account redshift distortions, whose effects on the 
power spectrum are at least as great as those of nonlinear- 
ity. There is no point in using these methods as they stand, 
without also taking into account redshift distortions. How- 
ever, given that a full-blown procedure including redshift 
distortions may well be based in part on the methods de- 
scribed, it seems worthwhile to lay out the steps required to 
implement them. 

9.1 Gourmet 

This version of the recipe is conceptually the simplest, but 
it takes the most computing power (a supercomputer would 
be handy). The procedure is a direct implementation of the 
FKP approximation to the minimum variance estimator of 
prewhitened power and the associated Fisher matrix, as de- 
scribed in §^. 

Naturally, if one were going to the trouble of using the 
gourmet recipe, then one would want to use the best possible 
model of the 3-point and 4-point correlation functions, not 
just the hierarchical model with constant amplitudes. 

Steps 1 and 2 below require knowledge of the selection 
function of a survey, but no actual data. Steps 3-5 require 
actual data from a galaxy survey. 

Step 1. Compute the FKP approximation to the asym- 
metric Fisher matrix E"^ of t he prewhitene d no n linea r 
power spec trum , as described in §S.l, equations (|l2C|)-(122). 



Equation (12C) involves a 5-dimensional integral over pairs 
ij of volume elements separated by r^j = r in the sur- 
vey. If the selection function n(r) separates into the prod- 
uct of an angular mask and a radial selection function. 



then the 3-dimensional angular integrals can be done an- 
alytically (Hamilton 1993), leaving a double integral of 
'^^\ka,ki3;ni,nj) over the radial directions. 
Step 2. The covariance matrix (AX q AXg) of the prewhi- 
tened power is equal, equation (116), to the inverse E^^^^^ 
of the symmetrized Fisher matrix. Use this covariance ma- 
trix to construct decorrelation matrices W , as described in 
Paper 4, with the property that W{AXAX^)W^ is diag- 
onal in Fourier space (cf. Paper 4 eq. [20]). The diagonal 
elements of this diagonal matrix are the expected variances 
of the decorrelated band-powers _B — WX to be computed 
in Step 5. 

Step 3. C omp u te th e estimator Za as described in §8.2 
equations ( 131 )-( 134 ). 
Step 4. Transform into the prewhitened pow er Xa using 
equations (||), (0), and (0), as stated in § |7.5| . 
Step 5. Decorrelate the estimated prewhitened power spec- 
trum X into a set of uncorrelated band-powers 13 = WX, us- 
ing the decorrelation matrices W computed in Step 2. Bear 
in mind that, as usual in ML fitting, the error bars should 
of course be interpreted as being attached to the model, the 
prior band-powers B, rather than to the data, the estimated 
band-powers 13. 



9.2 Fine 

This method adopts the approximation made in §^ and §^ 
that the prewhitened reduced covariance matrix 58 takes 
the simplified form (^) . According to the results of this 
approximation to should be quite good. If it is, then the 
fine method should yield results close to the gourmet method 
of §3.1, at a considerable saving in computer time. 

Steps 1-5 below do not require any actual data; the 
steps can be used to determine in advance how well the pre- 
whitened power spectrum might be measured from a survey. 
Steps 6-9 require actual data from a galaxy survey. 
Step 1. Compute a table of FKP-weighted pair integrals 
R{r;fi) at many separations r and several FKP constants 
/i. Calculating the pair integrals R{r; fi) requires knowing 
the selection function n(r) of a galaxy survey, but does not 
require actual data. This pair integral, commonly denoted 
{RR), is commonly computed by Monte Carlo integration, 
but I find it faster, more accurate, and more convenient 
(since the program was already written) to compute the in- 
tegral directly, using the procedures described by Hamilton 
(1993). 

Step 2. Compute the prewhitened 4-point contribution 
AI = H^^^^K H~^^^ to the reduced covariance of the non- 
linear power spectrum. This involves adopting a prior power 
spectrum ^{k), and a model of the 4-point correlation func- 
tion riijki. For the hierarchical model, the covariance matri- 
ces K and H are given by equations (^5|) and (^^. Some 
numerical issues concerning the computation of the matrix 



M are discussed in §4.1 



Step 3. Compute the eigenfunctions 4>a and eigenvalues 
equation (^^, by diagonalizing the prewhitened 4-point ma- 
trix M. 

Step 4. Compute the asymmetric Fisher matrix E"''^ , equa- 
tion (^), of the prewhitened nonlinear power spectrum of 
the survey, in the representation of eigenfunctions (l>a of the 
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prewhitened 4-point matrix M. This is where the pair inte- 
gral R{r, n) com pute d in Step 1 is needed. Numerical issues 



are discussed in i6.2. 

Step 5. Same as Step 2 of the gourmet method: from the 
inverse E'^^fs) '^^ symmetrized Fisher matrix, construct 
decorrelation matrices W such that the covariance of the 
band-powers B — WX is diagonal. Decorrelation is the sub- 
ject of Paper 4. 

Step 6. Compute a t able of FKP-weighted pair densities 
5'(r;/i), equation (102), at many separations r and several 



FKP constants /i. Calculating S{r;^), commonly denoted 
(DD) — 2{DR) 4- (RR), requires actual data from a survey. 
Step 7. From S{r;jj,a), compute the estimate Z'^ , equa- 
tion (105), in the representation of eigenfunctions (jia, as 
descri 



Ded in 



§0- 



step 8. Same as Step 4 of the gourmet method: transform 
Z°' into the prewhitened power Xc, using equations (p9|), 
([t^, and ([70|), as stated in § 7.5. The transformations may be 
done in whatever representation proves most convenient or 
numerically reliable. Ultimately, one wants the prewhitened 
power spectrum X{k) in Fourier space. 
Step 9. Same as Step 5 of the gourmet method: decorrelate 
the prewhitened power spectrum X into a set of uncorrelated 
band-powers B — WX, using the decorrelation matrices W 
computed in Step 5 above. 



9.3 Fastfood 

For some purposes a simplified, approximate version of the 
procedure in may be considered adequate. 

The basic simplifying approximation here is that the co- 
variance {AX{ka)AX{kf})) of the prewhitened power spec- 
trum may be considered to be diagonal in Fourier space 
without further refinement. The procedure then becomes the 
same as the FKP procedure for Gaussian fluctuations, with 
the differences that (a) it is the prewhitened power spectrum 
X{k), equation (p9[), rather than the power spectrum ^(fc) 
that is being estimated; and (b) the FKP constants ^(fc) 
in the FKP pair-weightings are modified from the Gaussian 
case where /j(fc) — S,{k). 

Figure fni shows the ratio ii{k)/£^{k) of the effective FKP 
constant /i(fcy to the nonlinear power spectrum ^(fc) for sev- 
eral different power spectra. The ratios plotted in Figure |l^ 
should be regarded as indicative rather than definitive, be- 
cause they depend on the validity of the hierarchical model 
with constant amplitudes Rb = —Ra (as required to satisfy 
the Schwarz inequality, §3.2), which as discussed in §^!^ is 
certainly wrong at some level. 

The effective FKP constant fj,{k) shown in Figure |ll| 
is not the same thing as the eigenvalue iia of Af at the 
nom inal wavenumber ka shown in Figure m As discussed 



4 , the correspondence between eigenvalue /la and nominal 



wavenumber ka is precise only for Gaussian fiuctuations. 

The effective FKP constant yu(fc) in Figure |ll] was cal- 
culated by going through Steps 2-5 of the recipe in j |9.2| for 
the case of a perfect, noiseless {fi —> <x) survey of (large) 
volume V. In this case the Fisher matrix Eap of the prewhi- 
tened power, equation (p3), reduces to 



V _i 



3.0 
2.5 
2.0 

mm) 

1.5 
1.0 




.0 L-L_L 



J LJ I LJ I I I L 
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Wavenumber k (/jMpc^') 

Figure 11. Ratio fj.{k)/^{k) of the effective FKP constant fi{k) to 
the nonlinear power spectrum ^(fe), as a function of wavenumber 
k, for several different power spectra. The ratios should be re- 
garded as indicative rather than definitive, because they depend 
on the validity of the hierarchical model (see text). The numbered 
curves are for power law power spectra with correlation functions 
^(r) = (r/5 h~^yipc) the number label being the index 7. The 
curve labelled P97 is for the Q^n = 0.3 power spectrum derived 
from observations by Peacock (1997), while that labelled ACDM 
is for the ACDM power spectrum of Eisenstein & Hu (1998). 



whose inverse gives the covariance of the prewhitened power 
X 



(AXcAXp) = - Map 



(137) 



Decorrelating this covariance, equation (137), in Fourier 
space, as described by in Paper 4, yields band-powers B{k) 
whose covariance is by construction diagonal. The diagonal 
values of the diagonal covariance matrix of the band-powers 
can be taken to define the effective FKP constants fi{k) 

2ll{kaf 



{AB{ka)AB{kp)) 



V 



(138) 



With the effective FKP constants fi{k) taken as given 
by Figure |l^, the shortcut recipe is then as follows. 
Step 1. Compute the effective spatial volume V{k) of the 
survey for modes at wavenumber k 

V(k) = 2,ikfRP,,m- j + ■ (139) 

Step 2. The Fisher matrix (^^ of the prewhitened power 
reduces to 

1 r°° 

E{ka,kfi) = - J jo{kar)jo{kf,r)R[r;fi{ka)]47vr'dr. (140) 

If the prewhitened power is averaged over sufficiently broad 
s hell s in fc-space, then, by arguments similar to those in 
jj6.3l the Fisher matrix is approximately diagonal (compare 



E{ka,k0) ^ {2TVy5:iD{ka 



V{ka) 
2^{ka) 



(141) 



(136) 



where V{ka) is the effective volume given by equation ( 139 ) 
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For the approximation (141) to be valid, the shells in fc-space 



must be broad not only compared to the inverse scale of the 
survey (as in the Gaussian case), but also compared to the 
width of the 4-point matrix M plotted in Figures ^ and ^ 
In the large k, hierarchical limit, the width of the matrix 
M in fc-space is comparable to an inverse correlation length, 
Afc ~ 7v/ro. 

Step 3. The covariance (AX(fca) AX(fe^)) of the pre- 
whitened power equals the inve rse of the Fisher matrix 
E{ka,ki3) given by equation (141) 



{AX{k^)AX{k0)) ^ {2nfS-iD{ko. - fc/j) 



V(fc„) 



(142) 



Define prewhitened band-powers B{k) to be the prewhitened 
power spectrum X{k) averaged over broad (as in Step 2) 
shells of volume Vk about k 



B{k) = 



i{k) dVfc 



(143) 



where dVfe = Ank^ Ak / (2-k)''^ . The variance of the shell- 
averaged prewhitened band-powers is 

2At(fc«)^ 



{AB{kY 



(144) 



V{k) Vk 

which is 2 ^{kaY divided by the effective phase volum e, th e 
product of the effective spatial volume V(k), equation (139), 
with the Fourier volume 14 of the shell in fc-space. 



Step 4. Same as Step 6 of §3.2: compute FKP- weighted pair 
densities S{r;fi), equation ( 102 ). 
Step 5. Compute the estimator Z(k) 



Mkr)S[r-Mk^ , 

' [i + C(^)]'/' 



(145) 



which may be compared to equation (104). 
Step 6. Same as^tep 8 of §9.2: transform Z{k) to X{k) 



usmg equations (§§), (0), and (0). The estimator Y = 
E-^Z is 



(146) 



Step 7. Form prewhitented band powers B{k) by aver- 
aging Xf fc) over sufficiently broad shells in fe-space, equa- 
tion (Fij). 



10 CONCLUSIONS 

The main finding of this paper is that the prewhitened non- 
linear power spectrum Xa defined by equation ( ]6£| ) has sur- 
prisingly sweet properties. 

Firstly, the covariance of the prewhitened nonlinear 
power is substantially narrower in Fourier space than the co- 
variance of the nonlinear power spectrum itself. Figures 

Secondly, in the FKP approximation, the 4-point and 
3-point contributions M and L to the covariance of prewhi- 
tened power are almost simultaneously diagonal (the 2-point 
contribution is by construction the unit matrix, so is auto- 
matically diagonal). Figures ^ and ^ Thus the eigenmodes 
of the covariance of prewhitened nonlinear power form a set 
of almost uncorrelated modes somewhat analogous to the 
Fourier modes of power in the Gaussian case. 



Thirdly, the eigenvalues Ha and A^, as defined by equa- 
tions (|6^) and (|68|), of the 4-point and 3-point prewhitened 
matrices AI and L are almost equal, fia ~ Aq, Figure ^ 
which is similar to the Gaussian case where jj,{k) = A(fc) = 

C(fc). 

The second and third points above together make it pos- 
sible to construct a near-minimum variance estimator, §^ 
and Fisher matrix, §^ of the prewhitened nonlinear power 
spectrum similar to the FKP estimator and Fisher matrix 
of the linear power spectrum in the Gaussian case. 

Fourthly, all the above properties hold for all power 
spectra tested, including power law nonlinear power spec- 
tra ^(fc) cx fc" with indices — 2 < n < over the full range 
allowed by the hierarchical model, and including realistic 
power spectra, such as the observationally derived power 
spectrum of Peacock (1997), and an observationally concor- 
dant ACDM model of Eisenstein & Hu (1998), nonlinearly 
evolved according to the Peacock & Dodds (1996) formula. 

Fifthly, in the realistic cases of the Peacock (1997) and 
Eisenstein & Hu (1998) power spectra, the prewhitened non- 
linear power spectrum X{k) appears to be curiously close to 
the linear power spectrum ^^(fc), Figure |l^. 

This having been said, it should be emphasized that 
the above properties are all premised on the hierarchical 
model with constant hierarchical amplitudes, §p.2|, which as 
discussed in 



4.2 



and by Scoccimarro et al. (1999 §3.3) is 
certainly wrong at some level. Clearly it will be important to 
test how well these results stand up in A'^-body simulations. 

In the meantime, the results of this paper raise ques- 
tions. Is there some physical reason underlying the seemingly 
unreasonably pretty properties of the prewhitened nonlinear 
power spectrum? In general, modes may be statistically un- 
correlated without being dynamically independent. But the 
fact that the covariance of the prewhitened power is nar- 
row for all power spectra is suggestive: do the eigenmodes 
of the covariance of prewhitened power somehow encode the 
information in the linear power spectrum that is ravelled 
by nonlinear evolution in the power spectrum itself? And 
is there somehow a connection to the mapping between lin- 
ear and nonlinear power spectra found by Hamilton et al. 
(1991)? 

I conclude with a repeat of the warning that this pa- 
per has ignored redshift distortions, light-to-mass bias, and 
evolution, and it has assumed that the only sources of vari- 
ance are cosmic variance and shot-noise variance arising 
from Poisson sampling of galaxies. In real galaxy surveys, 
all these problems must be grappled with. 
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APPENDIX A: JUSTIFICATION OF 
EQUATION (??) 



Equation ( |4l[ ) is the FKP approximation, expressed in con- 
cise mathematical form. This appendix offers further details 
justifying this equation. 

The pair-pair covariance matrix <tijki can be regarded 
as an operator that acts on pair-functions ^/fcj. It is help- 
ful to think of "i/ki as a 2-particle wavefunction (symmetric 
under pair exchange k ^ I), and Cijfci as a Hermitian op- 
erator that acts on the space of such wavefunctions. The 
pair-wavefunctions 'i!ki of interest in the present case have 
translation and rotation symmetry, which means that they 
have zero total momentum and zero total angular momen- 
tum. In the Fourier representation such wavefunctions 
can be expressed in the form 



*(fefc, ki) = {27vf S-ioikk + ki) ij{kk) 



(147) 



where ^{kk) is a function of th e scalar kk = jfefej. In a general 
representation, equation (147) is 



(148) 



where D^, is the operator introduced in §2.4, equations ( |28[ ) 
and (^. 

In the FKP approximation, the selection functions fii 
and fij upon which the pair-covariance ^ijki depends, equa- 
tions ( psj ) and (^5[), are taken to be locally constant, so that 
Cijfci also has translation and rotation symmetry, i.e. it com- 
mutes with the operators of total momentum and total an- 
gular momentum. Thus in the FKP approximation, Cijfci 
acting on a wavefunction D^iipa with zero momentum and 



angular momentum yields another wavefunction D^-xfi with 
zero momentum and angular momentum 



= D: 



Xl3 



(149) 



Take Tpa in equation (149) to be the elements o f a c omplete 
orthonormal basis of functions. Then equation (149) implies 
that 



D 



■I3a 



(150) 



for some matrix C^sq. Equation (150) is the desired equa- 
tion (|^) that was to be just ified (at least if the indices on 
£a/3 are swapped in eq. [150|, which is fine because is 
symmetric, as pr oven below). The wavefunctions tpa and xp 



in equation (149) are related by 



2/5^a = X0 



(151) 



A wavefunction of the form Dfjtpa is unnormalized - 
that is, ^lj°'D'JD^.'tpi3 diver ges - as is usual in quantum me- 
chanics for a wavefunction that has definite momentum, and 
must therefore be defined over infinite space. The divergence 
can be tamed by regarding the wavefunction as being defined 
instead over an extremely large but finite volume V. Then 



D'W'^ = 1 



a/3 



V 



(152) 



which is most easily proven from the real space represen- 



tation of Dfj, equation (pq). Equation (1152) s hould be in- 



terpreted with due care. For example, equation (152) should 
not be substituted into equation ( |45[ ) for the Fisher matrix in 
the FKP approximation, because the matrix {ni,nj) 
on the right hand side of equation ( p5| ) varies with positions 
i and j. 



Operating on equation (150) with D'^ implies, from 
equation (152), that (fii and fij here are being regarded for- 
mally as fixed constants in the huge volume V) 



13 a 



(153) 



It is evident from this equation that the pair exchange sym- 
metry ij ^ kl of Cijfcj implies that <tc,p is similarly sym- 
metric 



■pa 



(154) 



Equation (153) shows that, modulo the normalization fac- 



tor V , the reduced matrix ffa^ can be regarded as the ma- 
trix elements of the operator Cijfc; restricted to the class 
of wavefunctions that have zero total momentum and zero 
total angular momentum. 



APPENDIX B: FFTLOG 
B.l Introduction 

FFTLog computes the fast Fourier or Hankel {= Fourier- 
Bessel) transform of a periodic sequence of logarithmically 
spaced points. FFTLog can be regarded as a natural ana- 
logue to the standard FFT, in the sense that, just as the 
normal FFT gives the exact (to machine precision) Fourier 
transform of a linearly spaced periodic sequence, so also 
FFTLog gives the exact Fourier or Hankel transform, of arbi- 
trary order /i, of a logarithmically spaced periodic sequence. 
FFTLog shares with the normal FFT the problems of ring- 
ing (response to sudden steps) and aliasing (periodic folding 
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of frequencies), but under appropriate circumstances FFT- 
Log may approximate the results of a continuous Fourier or 
Hankel transform. 

The FFTLog algorithm was originally proposed by Tal- 
man (1978). However, it seems worthwhile here to present 
the algorithm in some detail. 

The FFTLog code may be downloaded from http://casa 
.colorado.edu/~ajsh/FFTLog/ . 

Consider the continuous Hankel (= Fourier-Bessel) 
transform pair 



a{r + R) = a{r) 



(162) 



a{k) — / a{r) (kr)'' Jfj,{kr) k dr 




poo 

a{r) — / a(fc) (fcr)~' Jp(fcr) r dfc . 
If the substitution 

a{r) = A{r) r^" and a(/c) = A{k) k" 



(155) 
(156) 

(157) 



is made, then the Hankel transform pair (|155|), (156), be- 
comes equivalent to the transform pair 



A{k) = / A{r)J^,{kr)kdr 



A{r) 



A{k)Jf,(kr)rdk 



(158) 



(159) 



Although the Hankel transform (155) with a power law bias 
(fer)*'' is thus equivalent in t he continuous case to the unbi- 
ased Hankel transform (158), the transforms are different 
when they are discretized and made periodic; for if a{r) 
is periodic, then A{r) — a{r) r'' is n ot p eriodic . FF TLog 
evaluates discrete Hankel transforms (155) and (156) with 
arbitrary power law bias. 

Fourier sine and cosine transforms can be regarded as 
special cases of Hankel transforms with /x = ±1/2, since 



-^1/2(2;) = (2/-KX) sm{x) 
J-i/2{x) ~ cos{x) 



(160) 
(161) 



As first noted by Siegman (1977), if the product kr 
in the Hankel transform is written as gin'^+i"'"^ then the 
transform becomes a convolution integral in the integration 
variable Inr or Infc. Convolution is equivalent to multiplica- 
tion in the corresponding Fourier transform space. Thus the 
Hankel transform can be computed numerically by the al- 
gorithm: FFT multiply by a function FFT back. This 
is the idea behind a number of Fast Hankel Transform algo- 
rithms (Candel 1981; Anderson 1982; Hansen 1985; Fanning 
1996) including FFTLog (Talman 1978). 

An advantage of FFTLog, emphasized by Talman 
(1978), is that the order jj, of the Bessel function may be 
any arbitrary real number. In particular, FFTLog works for 
1/2- integral /i, so includes the cases of Fourier sine and co- 
sine transforms, and spherical Hankel transforms involving 
the spherical Bessel functions jx{x) = {-k /2xY^'^ Jxj^i/2{x). 

B.2 Normal discrete Fourier transform 

First, recall the essential properties of the standard dis- 
crete Fourier transform of a periodic sequence of linearly 
spaced points. Suppose that a{r) is a continuous, in general 
complex-valued, function that is periodic with period R, 



Without loss of generality, take the fundamental interval to 
be [— _R/2, _R/2], centred at zero. Since a(r) is periodic, its 
continuous Fourier transform contains only discrete Fourier 
modes g^Timr/fl ^--^^Yi integral wavenumbers m. Suppose fur- 
ther that the function a{r) is 'smooth' in the specific sense 
that it is some linear combination only of the A'^ lowest fre- 
quency Fourier modes, m = 0, ±1, ±[A/2], where [A/2] 
denotes the largest integer greater than or equal to A/2, 



a(r) = Cm e 



(163) 



the outermost Fourier coefficients being equal, c_jv/2 ~ 
Cjv/2, in the case of even A. The primed sum in equa- 
tion (163) signifies a sum over integral m from —[A/2] to 
[A/2], with the proviso that for even A the outermost ele- 
ments of the sum receive only half weight: 



[JV/2] 

E ^ 

= -[JV/2] 



(164) 



with w„ = 1 except that iu_jv/2 = Wjv/2 = 1/2 if A is even. 

The sampling theorem (e.g. Press et al. 1986 §12. 1) as - 
serts that, given a function a(r) satisfying equation (163), 
the Fourier coefficients c™, can be expressed in terms of the 
values a„ = a(r„) of the function a{r) at the A discrete 
points r„ = nR/N for n = 0, ±1, ±[A/2]. For even A, 



the periodicity of a{r) ensures that a_ 



N/2 



Ojv/2- Specifi- 



cally, the sampling theorem asserts that the Fourier coeffi- 
cients in the expansion (163) satisfy 



J- -2-iTimn/N 

n2^ """^ 



the discrete points a„ themselves satisfying 



27riTnn / N 



(165) 



(166) 



in accordance w ith e quati on ( 163). 

Equations (165) and (166) constitute a discrete Fourier 
transform pair relating two periodic, linearly spaced se- 
quences a„ and Cm of length A. The standard FFT evaluates 
the discrete Fourier transform exactly (that is, to machine 
precision). 



B.3 Discrete Hankel transform 

Now suppose that the function a(r), instead of being pe- 
riodic in ordinary space r, is periodic in logarithmic space 
Inr, with logarithmic period L, 



a(re ) = a(r) 



(167) 



Take the fundamental interval to be [\nro~ L/2, lnro + L/2], 
centred at Inro. As in §B.2, the periodicity of a{r) im- 
plies that its Fourier transform with respect to In r contains 

2Trim}n{r / vq) / L 



only discrete Fourier modes e "^''/''oj/ wit h in tegral 
wavenumbers m. Suppose further, as in §B.2 eq. (163), that 
a(r) contains only the A lowest frequency Fourier modes 



a(r) = 



27rim ln{r / vq) / L 



(168) 
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with c 

serts that the Fourier coefficients Cm are given by 
1 



N/2 = C]v/2 for even A^. The samphng theorem as- 

(169) 



E -27rimn/JV 



where a„ = a(r„) are the values of the function a(r) at the 
N discrete points r„ = roe"^^^ for n = 0, ±1, ±[A^/2], 



E27rimn/JV 
Cm e 



(170) 



The continuous Hankel tra nsfo rm a{k), equation (155) 
of a function a{r) of the form (16S) is 



a{k) =^C^ g2™ln(r/ro)/L (^.^ ) " (fcr ) fc dr . (171) 

Jo 



The integrals on the right hand side of equation (171) can 
be done analytically, in terms of 

f/.(^)-Pv.Wd.^2^I||j±l±fj|| (172) 

wher e r (z) is the usual Gamma-function. Thus equa- 
tion ([l7l|) reduces to 



-27rimln(fc/fco)/L 



a{k) = ^ Cm? 

m 

where is 



9 + 



27rim\ 



-J 



(173) 



(174) 



Notice that = u-,„ , which ensures that a(fc) is real if a(r) 
is real. Equation ( |l73| ) gives the (exact) continuous Hankel 
transform a{k) of a function a{r) of the form (163). Like a{r), 



the Hankel transform a{k) is periodic in logarithmic space 
in k, with period L. The fundamental interval is [In fco — 1//2, 
in ko+L/2\, centred at In fco, which may be chosen arbitrarily 
(but see §B.5 below). 

The sampling theorem requires that u_jv/2 = un /q. fo r 
even A'^, which is not necessarily satisfied by equation ( 174 ) . 
However, at the discrete points fcn = koe"'"^^ considered 
by the sampling theorem, the contributions at m = ± N/2 
to the sum on the right hand side of equation (173) are 
( — )"e/v/2(n/v/9-|-ii t,/o)/ 2, whose imaginary part cancels out. 
Thus the equality (173) remains true at the discrete points 
fc„ if u±]v/2 are replaced by their real parts, 



U±]V/2 



Re 



Un/2 



(175) 



With the replacement (175), the sampling theorem asserts 
that the coefficients c-mUm in the sum (173) are determined 



by the values a„ = d{k„) of the Hankel transform at the 
discrete points k„ = fcoe"^/^ for n = 0, ±1, ±[N/2] 



Cm Um 



27vimn/ N 



\ ^ —27vimn/N 

0>7i — / CjnUm 6 



(176) 
(177) 



Putting together equations (|l69|), ( |l70| ), ( |l7q ) and ( |l77| ) 
yields the discrete Hankel transform pair 



(178) 



a™ = ^ a„i;„+„(/i,q) 



(179) 



in which the forward discrete Hankel mode v^{fi,q) is 
the di scret e Fou rier transform of Umiyi'jq) given by equa- 
tions dlTi) and dlTi), 



1 v-^' 

(M,g) = JJZ^ Um{fJ-,q)e 



27rimn / N 



(180) 



while the inverse discrete Hankel mode v„ {fi, q) is the dis- 
crete Fourier transform of the reciprocal l/u_m(pi, g), 



v-a (m, i) = 4; ]— 



-27rimn / N 



1) 



(181) 



The Hankel transform matrices v^^^{fi, q) and v^^„{fj,, q) 
are mutually inverse 



^ + ; (m, Q) «! + n (M, q) = Sr, 



(182) 



where Smn denotes the Kronecker delta. The forward and 
inverse Hankel modes have the interesting property of being 
self-similar; that is, Hankel modes v^^^(fi, q) [or v^+nif^j l)] 
with different indices m consist of the same periodic se- 
quence Vn{n, q) [or Vnin, q)] cyclically shifted by m notches. 

FFTLog evaluates the forward and i nv erse discrete Han- 
kel transforms given by equations (175), (17!:), exactly (to 
machine precision). 

The reciprocal l/it_m(A*, ?) in eq uati on ( flSlI) is equal 
to UmitJ-, ~q), according to equations (172) and (|174|), 



1 



.(M,-q) {ni^N/2) 



(183) 



except in t he c ase m — ±N/2 for even A'', w hen the re- 
placement (175) generally invalidates equation (183). How- 
ever, in t he sp ecial case where u±jv/2 are already real, then 
equation (175) leaves ii±jv/2 unchanged, and equation (18J) 
remains valid also at m = ±N/2. This special case is of 
particular interest, and is discussed further in §B.5 below. 

In the continuous case, the inverse Hankel transform is 
equal to the forward transform with q 



-g, equations (15 



and (156). In the discrete case this remains true for odd A'^, 
but it is not generally true for even A'^ (the usual choice) 

except in the important special case discussed in §B.5. 

In the general discrete case (i.e. if the condition [IS4] 
in §B.5 is not satis fied) , the inverse discrete Hankel mode 
v~{ij,,q), equation (181), d iffer s from the forward Hankel 
mode vi{jj,, —q), equation (180), only for even A' and only 



in the coefficient of the highest frequency Fourier compo- 
nent, l/u-m{fJ',q) versus UmifJ., —q) for m = ±N/2. To the 
extent that the highest frequency Fourier coefficient c±jv/2 
of a sequence a„ is small, the difference between its inverse 
discrete Hankel transform and its forward transform with 
q —q should be small. 

It is possible for the inverse discrete Hankel transform to 
be singular, if w±]v/2 is purely imaginary, so that its real part 
vanishes, making {/j,, q) singular. As discussed in §B.5, this 
singularity can be avo ided by choosing a low-ringing value 
of fcoro, equation ( |l8(j ) . 

The forward (inverse) discrete Hankel transforms are 
also singular at special values of /x and q, namely where 
+ 1 -I- q (or /i + 1 — g in the inverse case) vanishes, because 
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uo{fJ.,q) — U^{q) is singular at these points. This singular- 
ity reflects a real singularity in the corresponding continu- 
ous Hankel transform (unlike the singularity of the previous 
paragraph, which is an avoidable artefact of discreteness). 
The singularity in uq leads to an additive infinite constant 
in the discrete Hankel transform. In physical problems this 
additive infinite constant may somehow cancel out (for ex- 
ample, in the difference between two Hankel transforms). 
FFTLog's strategy in these singular cases is to evaluate the 
discrete Hankel transform with the infinite constant set to 
zero, and to issue a warning. 

B.4 FFTLog algorithm 

The FFTLog algorithm for taking the discrete Hankel trans- 
form, equation (17^), of a sequence a„ of TV logarithmically 
spaced points is: 

• FFT a„ to obtain the Fourier coefficients Cm, equa- 
tion (pig|); 



• multiply by Um given by equations (174) and (175) to 
obtain CmW^; 

• FFT CmUm back to obtain the discrete Hankel trans- 



form dn, equation (17'i 



A variant of the algorithm is to sandwich the above op- 
erations with power law biasing and unbiasing operations. 
For example, one way to take the unbiased continuo us H an- 
kel transform A{k) of a function A{r), equation (155), is 
to bias A{r) and A{k) with power laws, equation (157), and 
take a biased Hankel transform, equation (155). The discrete 
equivalent of this is: 

• Bias An with a power law to obtain a„ = Anr^'', equa- 
tion (157); 

• FFT a,i to obtain the Fourier coefficients Cm, equa- 
tion (pis|); 

• multiply by Um given by equations (174) and (175) to 
obtain CmUm', 

• FFT CmUm ba ck t o obtain the discrete Hankel trans- 
form a„, equation (177); 

• Unbias a„ with a power law to obtain An = a„kn'', 



equation (157) 



Although in the continuous limit the result would be iden- 
tical to an unbiased Hankel transform, in the discrete case 
the result differs. With a simple unbiased discrete Hankel 
transform, it is the sequence An that is taken to be peri- 
odic, whereas in the algorithm above it is not An but rather 
On that is periodic. 

The inverse discrete Hankel transform is accomplished 
by the same series of steps, except that Cm is divided instead 
of multiplied by Um- 

The FFTLog code is built on top of the NCAR suite of 
FFT routines (Swarztrauber 1979), and a modified version 
of an implementation of the complex Gamma-function from 
the gamerf package by Ooura (1996). 

FFTLog includes driver routines for the specific cases 
of the Fourier sine and cosine transforms. 

B.5 Low-ringing condition on fco^'o 

The central values Inro and Infco of the periodic intervals 
in Inr and Infc may be chosen arbitrarily. However, ring- 



ing of the discrete Hankel transform may be reduced, for 
either even or odd A*', if the product koro is chosen in such 
a way that the boundary points of the sequence Um, equa- 
tion (174), are equal 



-N/2 — Wjv/2 



(184) 



Recall that the general procedure, for e ven N, was to re place 
u±N/2 by their real part, equation (171: ). The condition ( 18'^ ) 
requires that k±jv/2 are already real. The condition (184) 
reduces ringing because it makes the periodic sequence Um 
fold smoothly across the period boundary at m = ±N /2. 

In addition to reducing ringing, the condition (184) 
means that equation (183) remains true also at m = ±N/2, 
so is true for all m. In this case the inverse Hankel mode 
Vn{yi,q)-, equation (181), is equal to the forward Hankel 



mode u„ (/i, — g) with q of the opposite sign 



{^^., q) = (m, -q) = -^y^ 



Um(^, -q) e 



-2'K\mn / N 



(185) 



In other words, if condition (184) is satisfied, then the 
inverse discrete Hankel transform equals the forward discrete 
Hankel transform with q —q. This is like the continuous 



Hankel transform, equations (155), ( |l5q ), where the inverse 
transform equals the forward transform with q —> ~q. 

The periodicity condition (184) on k±jv/2 translates, for 
real fj, and q, into a condition on fcoro 



ln(fco'"o) 



L_ 

N 



I^Arg Uf,{q + ^^^ +integer| (186) 



where Argz = Imlnz denotes the argument of a complex 
number, and integer is any integer. In other words, to reduce 
ringi ng, it may help to choose koro so as to satisfy the condi- 
tion (18f:). This is not too much of a restriction, since L/N 
is the logarithmic spacing betw een points {= one notch), so 
the low- ringing condition (186) allows koro to be chosen to 
lie within half a notch [= L/{2N)] of whatever number one 
chooses, for example within half a notch of koro = 1. 

FFTLog can be set to use automatically the low-ringing 
value of koro nearest to any input value of koro. 

How else does the choice of koro affect the Hankel trans- 
form? Increasing the value of \n{koro) by one notch L/N 
cyclic ally shifts the discrete Hankel transform a,i, equa- 
tion (177), by one notch to the left, a„ — > 5„_i. In other 



words, changing In(fcoro) by an integral number of notches 
shifts the origin of the transform, but leaves the transform 
otherwise unchanged, as might have been expected. 

In practice, since in most cases one is probably using the 
discrete Hankel transform as an approximation to the con- 
tinuous transform, one would probably want to use koro ~ 1 
(or 2, or vr, according to taste). 



B.6 Unitary Hankel transform 

The discrete Hankel transform with both low-ringing koro 
and no power law bias, g = 0, is of particular interest because 
it is unitary, like the Fourier transform. Indeed, being also 
real, the low-ringing unbiased Hankel transform is orthogo- 
nal, i.e. self-inverse, like the Fourier sine and cosine trans- 
forms. This is like the con ti nuou s unbiased (g = 0) Hankel 
transform, equations (155), ( [l5(^ ), which is self-inverse. 

The discrete Hankel modes Vm{^J^,Q) = Vm{^J^,Q) = 
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v^{fi,0) in the low-ringing unbiased (q — 0) case are pe- 
riodic, orthonorrnal, and self-similar, equation (182), 



I 



(187) 



Like any orthogonal transformation, the low-ringing un- 
biased {q = 0) Hankel transform commutes with the opera- 
tions of matrix multiplication, inversion, and diagonalization 
(for non-low-ringing or biased Hankel transforms, q ^ 0, the 
operations do not commute). That is, the Hankel transform 
of the product of two matrices is equal to the product of 
their Hankel transforms, and so on. 

All else being equal (which it may not be), given a choice 
between applying an unbiased {q — 0) or biased (q 0) 
Hank el tr ansform, and between a low-ringing fcoro, equa- 
tion (186), or otherwise, one would be inclined to choose 
the low-ringing unbiased transform, because of its orthogo- 
nality property. 



B.7 Example 

Figure |l2| shows the correlation function ^(r) computed by 
FFTLog for the nonlinear ACDM power spectrum of Eisen- 
stein & Hu (1998) shown in Figure ^ Two different resolu- 
tions are plotted on top of each other, a low resolution case 
with 96 points over the range r — 10~^ to 10'' h~^Mpc, and a 
high resolution case with 768 points over the range r = 10~® 
to 10*^ /i"^Mpc. Both cases used an unbiased {q = 0) trans- 
form and a low-ringing value of koro (actually the choice of 
kovo made little difference here). 

The low and high resolution correlation functions shown 
in Figure p_a agree well except near the edges r ~ 10"'^ 
and 10^ h~ Mpc; in particular, the low resolution corre- 
lation function tends to a positive constant « 10~^ at 
r — > 10^ /i~^Mpc, whereas the high resolution correlation 
function is negative and declining as a power law cx r"* at 
large r. The disagreement is caused by aliasing (see §B.8) of 
small and large separations in the low resolution case. Alias- 
ing is almost eliminated in the high resolution case because 
the ranee r over which the transform 

was computed is much broader than the range plotted. 

The bottom panel of Figure ^ shows the ratio 
CpFT/CpFTLog of the correlation function ^fft computed 
with a normal FFT (sine transform) with 1023 points over 
the range r — 0.125 to 128ft~^Mpc, to the (high resolu- 
tion) correlation function ^FFTLog computed with FFTLog. 
Even with 1023 points, the FFT'd correlation function rings 
noticeably, with an amplitude of about ±5 percent. 

In this particular instance, FFTLog outperforms the 
normal FFT on all counts: it is more accurate, with fewer 
points, over a larger range, and it shows no signs of ring- 
ing. This does not mean that FFTLog is always better than 
FFT. Rather, FFTLog is well matched to the problem at 
hand: the cosmological power spectrum extends over many 
orders of magnitude in wavenumber k, and varies smoothly 
in In k. 
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Figure 12. Correlation function ^(r) corresponding to the non- 
linear COBE-normalized ACDM power spectrum of Eisenstein &c 
Hu (1998), Figure The top panel shows the correlation func- 
tion computed with FFTLog at two different resolutions, plotted 
on top of each other: (a) with 96 points over the range r = 10~^ 
to 10'^ /i~^Mpc (low resolution), and (b) with 768 points over 
the range r = 10-^ to lO'^h-^Mpc (high resolution). The lines 
are dashed where the correlation function is negative, at separa- 
tions r > 119h~^Mpc. The low and high resolution curves are 
almost indistinguishable except at r ^ 200/i~'^Mpc: the high res- 
olution curve is the one that declines as a power law ~ r~* at 
large r. The straight dashed line shows the canonical power law 
{r /5 h~^Mpc)~^'^ for reference. The middle panel shows the ra- 
tio 5iow/5high of the low to high resolution correlation functions. 
The bottom panel shows the ratio ^FFT/^FFTLog of the correla- 
tion function ^fft computed with a normal FFT (sine transform) 
with 1023 points over the range r = 0.125 to 128 h~^Mpc, to the 
high resolution correlation function ^FFTLog computed with FFT- 
Log. The FFT'd correlation function ^fft rings at the ±5 percent 
level. 



B.8 Ringing and aliasing 

FFTLog suffers from the same problems of ringing (response 
to sudden steps) and aliasing (periodic folding of frequen- 
cies) as the normal FFT. 

Usually one is interested in the discrete Fourier or Han- 
kel transform not for its own sake, but rather as an approx- 
imation to the continuous transform. The usual procedure 
would be to apply the discrete transform to a finite segment 
of the function a(r) to be transformed. For FFTLog, the 
procedure can be regarded as involving two steps: truncat- 
ing the function to a finite logarithmic interval, which causes 
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ringing of the transform; followed by periodic replication of 
the function in logarithmic space, which causes aliasing. 

Figure ^ illustrates these steps for the unbiased {q = 0) 
Hankel transform, equation (155), of order fi = —1/2 of a 
function that is Gaussian in the log 



Function a(r) Transform d(k) 

r*io-'iO"^iO"' I 10 10^ 10-' 10* io-'io-3io-2io-> i lo lo- lo' lo* 



i(r) = exp[-(lnr-)V2] 



(188) 



Truncation of the function a{r) leads to ringing of its 
transform a{k) at high frequencies k, as seen in the middle 
right panel of Figure The oscillations at large k are actu- 
ally uniformly spaced in k, but appear bunched up because 
of the logarithmic plotting. 

Periodic replication means taking a sum of copies 
shifted by integral periods. From the definition (155) of the 
continuous Hankel transform, it can be seen that period- 
ically replicating a function a(r) in logarithmic space Inr 
and then taking its continuous Hankel transform is equiva- 
lent to Hankel transforming the function a(r) and then peri- 
odically replicating the transform a{k) in Infc. But truncat- 
ing a function does not truncate its transform. So whereas 
a truncated, periodically replicated function a(r) contains 
contributions from only one period at each point r, the pe- 
riodically replicated transform contains overlapping contri- 
butions from many periods at each point k. This is aliasing. 
In Figure |l^ aliasing is visible as an enhancement of the pe- 
riodically replicated transform d{k) on the high k side of the 
periodic interval. 

Ringing and aliasing can be reduced by taking suitable 
precautions. 

The ringing that results from taking the discrete trans- 
form of a finite segment of a function can be reduced by ar- 
ranging that the function folds smoothly from large to small 
scales. It may help to bias the function with a power law 
before transforming it, as in the second algorithm in §B.4. 
It may also help to use a low-ringing value of fco^o, §B.5. 

Aliasing can be reduced by enlarging the periodic inter- 
val. Aliasing can be eliminated (to machine precision) if the 
interval can be enlarged to the point where the transform 
a{k) goes sensibly to zero at the boundaries of the period. 
Note that it is not sufficient to enlarge the interval to the 
point where a(r) is sensibly zero at the period boundaries: 
what is important is that the transform a{k) goes to zero at 
the boundaries. 
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Figure 13. Illustrating the ringing and aliasing that occurs when 
the continuous Hankel transform of a function is approximated by 
the discrete Hankel transform of a finite segment of the function. 
Lines are dashed where values are negative. The function a(r) is 
shown to the left, and its corresponding Hankel transform a{k) to 
the right. The panels from top to bottom are: (top) the original 
function a(r) and its Hankel transform d{k); (middle) the trun- 
cated function a(r) and its Hankel transform a(fc), which rings 
at high frequencies k; and (bottom) the truncated, periodically 
replicated function a(r) and its corresponding periodically repli- 
cated Hankel transform a(fc), which is aliased. Vertical lines in 
the bottom panels demarcate periodic intervals. 
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